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1 Introduction 



The ultrashort laser pulses, i.e. the pulses with the durations ~10~ 10 -1CP 15 
sec, have a lot of the applications, which range ultrafast spectroscopy, tracing 
chemical reactions, precision processing of materials, optical networks and com- 
puting, nuclear fusion and X - ray lasing, ophthalmology and surgery (for review 
see ly]). The mechanisms of the ultrashort pulse generation are active or passive 
loss switching (so-called Q-switching, Part II) and locking of the longitudinal 
laser modes (Part III) due to the active (Part V) or passive ultrafast modulation 
resulting in the laser quasi-soliton formation. Such quasi-soliton is very similar 
to the well-known Schrodinger soliton, which runs in the optical networks (Part 
VI). As a matter of fact, the model describing active mode locking is based 
on the usual equation of the harmonic oscillator or its nonlinear modifications, 
while for the passive mode locking description we have primordially nonlinear 
Landau-Ginzburg equation (Part VII). This equation is the dissipative analog 
of the nonlinear Schrodinger equation and, as a result of the nonlinear dissipa- 
tion, there exist a lot of nonstationary regimes of the ultrashort pulse generation 
(Part VIII). This requires the generalization of the model, which leads to the 
numerical simulations on the basis of FORTRAN (or C) codes generated by 
Maple (Part IX). Simultaneously, the obtained numerical results are supported 
by the analytical modelling in the framework of the computer algebra approach. 
The last takes into account the main features of the nonlinear dissipation in the 
mode-locked laser, viz. the power- (Part VII) or energy-dependent (Part X) 
response of the loss to the generation field. In the latter case, there is the pos- 
sibility of the so-called self-induced transparency formation, which is described 
by the nonlinear Klein-Gordon equation (Part XI). 

Our considerations are based on the analytical or semi-analytical search of 
the steady-state soliton-like solutions of the laser dynamical equations and on 
the investigation of their stability. Also, the breezer-like solutions are consid- 
ered using the aberrationless approximation. The computer algebra analysis is 
supported by the numerical simulations on the basis of the Maple generated 
FORTRAN-code. We present the analysis of these topics by means of the pow- 
erful capacities of Maple 6 in the analytical and numerical computations. This 
worksheet contains some numerical blocks, which can take about of 12 Mb and 
18 min of computation (1 GHz Athlon). 



2 Nonstationary lasing: passive Q-switching 



2.1 Continuous- wave oscillation 

The basic principle of Q-switching is rather simple, but in the beginning let's 
consider the steady-state oscillation of laser. The near-steady-state laser con- 



taining an active medium and pumped by an external source of the energy 
(lamp, other laser or diode, for example) obeys the following coupled equations: 
>restart: 
>with(linalg): 

>eql := diff(Phi(t),t) = (alpha(t) - rho)*Phi(t);# field evolution 
>eq2 := diff(alpha(t),t) = sigma[p]*(a[m]-alpha(t))*P/(h*nu[p]) - 
alpha(t)*sigma[g]*Phi(t)/(h*nu[g]) - alpha(t)/T[r];# gain evolution for quasi- 
two-level active medium 



eql :=^*(t) = (a(t)-p)*(t) 
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Here <&(£) is the time-dependent field intensity, a(t) is the dimensionless gain 
coefficient, P is the time- independent (for simplicity sake) pump intensity, v p 
and v p are the frequencies of the pump and generation fields, respectively, a p 
and a g are the absorption and generation cross-sections, respectively, T r is the 
gain relaxation time, p is the linear loss inclusive the output loss of the laser 
cavity, and, at last, a m is the gain coefficient for the full population inversion 
in the active medium. The pump increases the gain coefficient (first term in 
eq2), that results in the laser field growth (first term in eql). But the latter 
causes the gain saturation (second term), which can result in the steady-state 
operation (so-called continuous- wave, or simply cw, oscillation): 

>rhs( subs({alpha(t)=alpha,Phi(t)=Phi},eql) ) = 0; 
>rhs( subs({alpha(t)=alpha,Phi(t)=Phi},eq2) ) = 0; 
>sol := solve({%,%%},{Phi,alpha}); 
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The second solution defines the cw intensity, which is the linear function of 
pump intensity: 

>expand( subs( sol[2],Phi ) ): 

>Phi = collect(%,P); 

>plot( subs( {lambda[g]=8e-5, lambda[p]=5.6e-5}, 
subs( {h=6.62e-34, sigma[g]= 3e-19, sigma[p]= le-19, rho=0.1, nu[g]=3el0/ 
lambda[g], nu[p]=3el0/lambda[p], T[r]=3e-6, a[m]=2.5},rhs(rho*%) ) ), P = 
4.9e4..1e5, axes=boxed, title='output intensity vs. pump, [W/cm 2 ]' );# lambda 
is the wavelength in [cm] 
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The pump corresponding to $=0 defines so-called generation threshold. Now 
let's consider the character of the steady-state points sol of our system {eql, 
eq2}. The Jacobian of the system {eql, eq2} is: 

>eq3 := subs( {Phi(t)=x,alpha(t)=y},rhs(eql) ):# x = Phi(t), y = alpha(t) 
>eq4 := subs( {Phi(t)=x,alpha(t)=y},rhs(eq2) ): 

>A := vector( [eq3, eq4] );# vector made from the right-hand side of 
system {eql, eq2} 

>B := jacobian(A, [x,y]); 
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For the cw-solution the eigenvalues of the perturbations are: 

>BB := eigenvalues(B): 

>BBB := subs( {x=subs( sol[2],Phi ),y=subs( sol[2],alpha )},BB[1] ): 

>BBBB := subs( {x=subs( sol[2],Phi ),y=subs( sol[2], alpha )},BB[2] ): 
>plot( [subs( 
{lambda[g]=8e-5, lambda[p]=5.6e-5} 
,subs( 

{h=6.62e-34, sigma[g]= 3e-19, sigma[p]= le-19, rho=0.1, nu[g]=3el0/lambda[g], 
nu[p]=3el0/lambda[p], T[r]=3e-6, a[m]=2.5} 
,BBB)), 
subs( 

{lambda[g]=8e-5, lambda[p]=5.6e-5} 
,subs( 

{h=6.62e-34, sigma[g]= 3e-19, sigma[p]= le-19, rho=0.1, nu[g]=3el0/lambda[g]. 
nu[p]=3el0/lambda[p], T[r]=3e-6, a[m]=2.5} 

,BBBB/le7 ) )], P=4.9e4..1e5, axes=boxed, title='perturbation eigenvalues vs. 
pump' );# second eigenvalue is divided by le7 
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So, cw oscillations is stable in our simple case because the eigenvalues are 
negative. For the zero field solution the perturbation eigenvalues are: 

>BB := eigenvalues(B): 

>BBB := subs( {x=subs( sol[l],Phi ),y=subs( sol[l],alpha )},BB[1] ): 

>BBBB := subs( {x=subs( sol[l],Phi ),y=subs( sol[l],alpha )},BB[2] ): 
>plot( [subs( 
{lambda[g]=8e-5, lambda[p]=5.6e-5} 
,subs( 

{h=6.62e-34, sigma[g]= 3e-19, sigma[p]= le-19, rho=0.1, nu[g]=3el0/lambda[g], 
nu[p]=3el0/lambda[p], T[r]=3e-6, a[m]=2.5} 
,BBB ) ), 
subs( 

{lambda[g]=8e-5, lambda[p]=5.6e-5} 
,subs( 

{h=6.62e-34, sigma[g]= 3e-19, sigma[p]= le-19, rho=0.1, nu[g]=3el0/lambda[g], 
nu[p]=3el0/lambda[p], T[r]=3e-6, a[m]=2.5} 

,BBBB/le7 ) )], P=4.9e4..1e5, axes=boxed, title='perturbation eigenvalues vs. 
pump' );# second eigenvalue is divided by le7 
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The existence of the positive eigenvalue suggests the instability of the zero- 
field steady-state solution. Hence there is the spontaneous generation of the cw 
oscillation above threshold in the model under consideration. 



2.2 Q-switching 

The situation changes radically due to insertion of the saturable absorber into 
laser cavity. In this case, in addition to the gain saturation, the loss saturation 
appears. This breaks the steady-state operation and produces the short laser 
pulses. 



As a result of the additional absorption, Q-factor of laser is com- 
paratively low (high threshold). This suppresses the generation. 
When $ is small, the gain increases in the absence of the gain 
saturation (see eq2 from the previous subsection). This causes 
the field growth. The last saturates the absorption and abruptly 
increases Q-factor. The absorption "switching off" leads to the 
explosive generation, when the most part of the energy, which is 
stored in the active medium during pumping process, converts 
into laser field. The increased field saturates the gain and this 
finishes the generation. 



As the reference for the model in question see, for example, ||. To formulate 
the quantitative model of the laser pulse formation let's use the next approx- 
imations: 1) the pulse width is much larger than the cavity period, and 2) is 
less than the relaxation time, 3) the pump action during the stage of the pulse 
generation is negligible. We shall use the quasi-two level schemes for the gain 
and loss media (the relaxation from the intermediate levels is fast). Also, the 
excited-state absorption in absorber will be taken into account. 

The system of equation describing the evolution of the photon density (j)(t) 
is 

>restart: 

>with(plots): 

>print('System of basic equations:'); 

>el := Diff(n[l](t),t) = 
-sigma[s]*c*phi(t)*n[l](t);# EVOLUTION OF ABSORPTION. The ground level 
population n[l](t) defines the absorption (quasi-two level scheme). The relax- 
ation is slow in the comparison with the pulse duration, phi(t) is the photon 
density 

>e2 := Diff(phi(t),t) = 
(phi(t)/t_cav)* 

(2*sigma[g]*x(t)*l - log(l/R) - L - 2*sigma[s]*n[l](t)*l[s]-2*sigma[esa]*(n[0]- 
n[l](t))*l[s]);# EVOLUTION OF PHOTON DENSITY. The field variation over 
cavity round-trip is small, sigma[g] is the gain cross-section, x(t) is the inver- 
sion in amplifier defining the gain coefficient, t_cav is the cavity period, 1 is the 
active medium length, R is the output coupler refractivity, L is the linear loss, 
l[s] is the absorber length 

>e3 := Diff(x(t),t) = 
-gamma*sigma[g]*c*phi(t)*x(t);# EVOLUTION OF GAIN, c is the light veloc- 



ity, gamma is the parameter of the inversion reduction (2 for the pure three-level 
scheme and 1 for the pure four-level scheme) 



System of basic equations : 



d 

el := — m (t) = -o s c 4>{t) m (t) 



o <p(t) (2a g x(t)l -ln(— ) - L- 2a s n 1 (t)l s - 2a esa (n -ni(t))l s ) 

e2 := ^ 4>{t) = ^ 

at t cav 



d_ 

di 



e3 := — x(t) = -7 <j g c <j)(t) x(t) 



Now we shall search the ground state population in the absorber as a function 
of the initial population inversion in amplifier: 

>Diff(nl(n),n) = 
subs({n[l](t)=n[l](n),x(t)=x},rhs(el)/rhs(e3));# devision of el by e3 

>lnt(l/y,y=n[0]..n[l]) = Int(zeta/z,z=x_i..x);# zeta = sigma[s]/gamma/sigma[g], 
x_i is the initial inversion in amplifier, n[0] is the concentration of active ions 
in absorber 

>solve(value(%),n[l]): 
>n[l] = expand(%); 

>print('Ground state population in absorber:'); 
>sol_l := n[l] = n[0]*(x/x_i)^; 
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nl(n) 


a s ni(n) 


dn 
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Ground state population in absorber 



sol 1 := ri\ = no ( p 

x i 



The similar manipulation allows to find the photon density as a function of 
inversion in amplifier: 

>A := Diff(phi(x),x) = 
simplify ( subs( n[l](t)=rhs(sol_l),subs( {x(t)=x},rhs(e2)/rhs(e3) ) ) );# divi- 
sion of e2 by e3 

>B := numer( rhs(A) ); 
>C := denom( rhs(A) ); 
>BB := subs( 
{op(4,B) = (x/x_i)«*log(l/T[0] 2 ), 
2*sigma[esa]*l[s]*n[0]=ln(l/T[s] 2 ), 

op(6,B)=-ln(l/T[s] 2 )*(x/x_i)^}, B );# T[0] is the initial transmission of the 
absorber (log(l/T[0] 2 )= 2*sigma[s]*l[s]*n[0], l[s] is the absorber length) 

>CC := 2*l_cav*gamma*sigma[g]*x;# l_cav=t_cav*c/2 is the cavity length, 
t_cav is the cavity period 



A:=-H-<f>(x) = 

ox 

I Op r£ 

-2cr„x/ + ln( — ) + L + 2cr s n ( : )£ l s + 2cr esa l s n - 2cr esa l s n ( : ) c 

K x _i x _i 

t_cav 7 o g ex 



B := -lugxl + ln(— ) + L + 2a s n ( : ) c l s + 2a esa l s no - ^o- esa l s n ( : ) c 

K xi xi 



C := t_cav 7<r s ex 



BB:=-2a g xl + Mh + L + (-^)( ln(-L) + ln(-L) - ln(-L) (— )< 

R x _i Tq T s T s x _% 



CC := 2 l_cavj a g x 



>print ('Evolution of the photon density:'); 
>e4 := diff(phi(x),x) = 



10 



subs(ln(l/(T[s] 2 ))=delta*ln(l/(T[0] 2 )),BB)/CC;#delta=sigma[esa]/sigma[s] = 
ln(T[s])/ln(T[0]) is the parameter defining the contribution of an excited-state 
absorption with cross-section sigma[esa], T[s] is the fully saturated transmission 
of absorber 



Evolution of the photon density 



e 4 := ^ m = 



-<51n(^)( — )< 
T 2 ' x_i> 



l_cav j o~ g x 



Hence the photon density is: 



>dsolve({e4,phi(x_0)=0},phi(x)): 
>simplify( subs(x_0=x_i,%) ): 

>print('This is the basic dependence for photon density:'); 

>sol_2 := collect(combine(%,ln),{log(l/T[0] 2 ),sigma[g],zeta}); 



This is the basic dependence for photon density 



, N 1 2xl-2lx i 

sol_2 := 4>{x) = : — - 

2 l_cavj 



1 5 (ln(a;) - ln(a;_i)) ln(— -%) 1 -ln(-j-) (ln(cc) - ln(x_i)) - L (ln(cc) - ln(x_i)) 



l_cav 7 



2 



l_cav 7 



( _ ( ^)C + 5( ^)C_ ( j + i )ln( J_ ) 

1 x _i x _% To 

2 l_cav 7 er g £ 
So, we have: 
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' cav 7 



(Eq. 1) 
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Now let's define the key Q-switching parameters: 

>subs( n[l](t)=n[0],rhs(e2)*t_cav/phi(t)) = 0;# Q-switching start, n[l](0) 
= n[0] 

>print('Solution for the initial inversion:'); 

sol_3 := x_i = subs( 2*sigma[s]*n[0]*l[s]=log(l/T[0] 2 ), solve(%,x(t)) ); 



2 a g x(t) I — ln( — ) — L — 2<7 s nol s = 
R 



Solution for the initial inversion : 



1 '4» + i + " 1 <^» 

sol 3 := x i = 



o g l 



So, the initial inversion defining the gain at Q-switching start is 



ln(i)+ln( i )+L 



2a g l 



(Eq. 2) 



>e5 := numer( rhs(e4) ) = 0;# definition of the pulse maximum 
> print ('The inversion at the pulse maximum:'); 
>sol_4 := x_t = solve(e5,x); 



e 5:=-2a s ^ + ln(i)+L+(^ln(^)+51n(-L)-51n(-L)(-^)C = 



The inversion at the pulse maximum : 



sol_4 '■= x _t = 

RootOf (2 e- z I x lCTg -ln(i)-L-e(- zc >ln(-i T )- I 51n( 1 )+51n( 1 )e<- z «>) 

e ~~ T ° T ° T ° xi 
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>e6 := numer( simplify( rhs(sol_2) ) ) = 0;# definition of the Q-switching 
finish 

>print('The inversion at Q-switching finish'); 
>sol_5 := x_f = solve(e6,x); 



e6 := -2 1 x cr g ( + 2 1 x _i cr g ( + ln(x) 5\n( 7 ^) < + ln(:r) ln( — ) ( + ln(x) L ( 



v i I )C + ln(x)ln(i 



- ]n(x_i)6ln(^) ( - ln(x_i) ln(i) C - ln(x_i) L ( + (— )< ln(-L) 

- *M^ff) (^-) c + SH^s) H^s) = o 

The inversion at Q — switching finish 
sol _5 :— x _f = 

RootOi{2e- z lx_ia g C-2lx_icr g (~51n{jl^)( _Z-ln(j ? )( _Z-L( _Z-e<-- z ° ln(^) 
To^' To z T a z x _l 

Additionally, we define the inversion at the pulse maximum, when £ tends 
to infinity: 

>subs( (x/x_i)^=0,e5 );# by virtue of x_i > x and zeta — > infty 
> print ('Inversion at the pulse maximum for large zeta'); 
>sol_6 := x_tO = solve(%,x); 



-2a g xl + ln(-) + L + 6ln(— 
-ft i n 



Inversion at the pulse maximum for large zeta 



ln(i) + L + <y]n( * ) 
sol 6:=x tO- ' ° 



2 <7 S / 



So, we have the expressions for Xi (initial inversion, sol_3), Xt (the inversion 
at pulse maximum, sol_4 and e5), Xt t o (the inversion at pulse maximum when 
£ — ► oo, sol_6), Xf (the final inversion, sol_5 and e6) and the photon density 
4> as function of inversion x (sol_2). 
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As an example, we consider the real situation of Yb/Er-glass laser with the 
crystalline Co:MALO saturable absorber. The obtained expressions allow to 
plot the typical dependencies for the pulse parameters: 

>print('Pulse energy:'); 

>fun_l := (h*nu*S)/(2*sigma[g]*gamma)*log(l/R)*log(x_i/x_f); 
>print('Pulse power:'); 

>fun_2 := (h*nu*S*l)/(t_cav*gamma)*log(l/R)*(x_i- x_t - x_tO* 
log(x_i/x_t) - (x_i-x_tO)*(l-(x_t/x_i) a )/a); 
>print('Pulse width:'); 

>fun_3 := simplify (fun_l/fun_ 2); 



Pulse energy : 






Pulse power : 



huSlln(-) 
K 



fun_2 



( ■ (x_i - x_tO) (1 - (^4) a A 

x i — x t — x W ln( ~ ) = 

X _t a 

V J 

t _ cav 7 



Pulse width 



fun_3 := — 

x i 
ln( — =-) t cav a 
1 x_f 



O X % X t X t 

a g l(-x_ia + x_ta + x_W ln(^^) a + x_i - x_i (^^-) a - x_tO + x_tO (^^-) a ) 



This numerical procedure plots the dependence of the output pulse energy 
on the reflectivity of the output mirror: 
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>En := proc(gam,L,T_0,l,wO,i) # definition of system's parameters 
>delta := 0.028: 

>sigma[g] := 7e-21:# the gain cross-section in [cm 2 ] 

>sigma[s] := 3.5e-19:# the absorption cross-section in Co:MALO 

> alpha := sigma[s]/(sigma[g]*gam): 
>h := 6.63e-34:# J*s 
>nu := evalf(3e8/1.535e-6):# lasing frequency for 1.54 micrometers 
>S := Pi*w0 2 /2:# area of Gaussian beam in amplifier, wO is the beam ra- 
dius 
>R := 0.5+0.5*i/100: 

>x_i := l/2*(ln(l/R) + L + ln(l/(T_0 2 )))/(l*sigma[g]): 
>eq := 
ln(x)*delta*ln(l/(T_0 2 ))*alpha - ln(x_i)*delta*ln(l/(T_0 2 ))*alpha - 2*1* x*sigma[g]* 
alpha + 2*l*x_i*sigma[g]*alpha + ln(x)*L* alpha + ln(x)*ln(l/R)*alpha - 
ln(x_i)*L*alpha - ln(x_i)*ln(l/R)*alpha + (x/x_i) Q *ln(l/(T_0 2 )) - delta* 
ln(l/(T_0 2 ))*(x/x_i) Q - ln(l/(T_0 2 )) + delta*ln(l/(T_0 2 )) = 0:# e6 

>sol_f := fsolve(eq, x, avoid={x=0}): 

>sol_En := evalf( l/2*h*nu*S*ln(l/R)*ln(x_i/sol_f)/(sigma[g]*gam) 
* le3 ):# [mJ] 

end: 

>print('The parameters are:'); 
>print('l) inversion reduction factor gamma'); 
>print('2) linear loss L'); 

>print('3) initial transmission of absorber T[0]'); 
>print('4) gain medium length 1 in cm'); 
>print('5) beam radius in amplifier wO in cm'); 



The parameters are : 



1 ) inversion reduction factor gamma 



2) linear loss L 



3) initial transmission of absorber T[0] 
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4) gain medium length I in cm 



5) beam radius in amplifier wO in cm 



For the comparison we use the experimental data (crosses in Figure): 

>points := 
{seq([0.5 + 0.5*k/100,En(1.9,0.04,0.886,4.9,0.065,k)],k=1..100)}: 

>points_exp := {[0.793, 10.5], [0.88,9], [0.916,5. 5]}: 

>plot(points,x=0.5..1, style=point, axes=boxed, symbol=circle, 
color=black, title='Pulse energy vs. R, [mJ]'): 

>plot(points_exp,x=0.5..1, style=point, symbol=cross, 
axes=boxed, color=red): 
>display({%,%%}); 
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Pulse energy vs. R, [mJ] 
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Similarly, for the output power we have: 

>Pow := proc(gam,L,T_0,l,l_cav,wO,i) # definition of system's parameters 
>delta := 0.028: 

>sigma[g] := 7e-21:# in [cm 2 ] 
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>sigma[s] := 3.5e-19:# Co:MALO 

> alpha := sigma[s]/(sigma[g]*gam): 
>h := 6.63e-34:# J*s 
>nu := evalf(3e8/1.354e-6):# frequency for 1.54 micrometers 
>S := Pi*w0 2 /2:# area for Gaussian beam 
>R := 0.5 + 0.5*i/100: 
>c := 3el0: 

>refractivity := 1.6:# refractivity coefficient for the active medium 
>t_cav := 2*(l_cav - l)/c + 2*(l*refractivity)/c: 

>x_i := l/2*(ln(l/R) + L + ln(l/(T_0 2 )))/(l*sigma[g]): 

>x_t0 := 
l/2*(ln(l/R) + L + delta*ln(l/(T_0 2 )))/(sigma[g]*l): 

>eq := 
-2*sigma[g]*x*l + ln(l/R) + L + (x/x_i) a *ln(l/(T_0 2 )) + delta*ln(l/(T_0 2 )) 
- delta*ln(l/(T_0 2 ))*(x/x_i) Q = 0:# e5 

>sol_t[i] := fsolve(eq, x, avoid={x=0}): 
>sol_Pow[i] := h*nu*S*l*ln(l/R)* 
(x_i-sol_t[i] - x_tO*ln(x_i/sol_t[i]) - (x_i-x_tO)*(l-(sol_t[i]/x_i) a ) /alpha)/ 
(t_cav*gam)/le3:# [kW] 

>end: 

>print('The parameters are:'); 
>print('l) inversion reduction factor gamma'); 
>print('2) linear loss L'); 

>print('3) initial transmission of absorber T[0]'); 
>print('4) gain medium length 1 in cm'); 
>print('5) cavity length l_cav in cm'); 
>print('6) beam radius in amplifier wO in cm'); 



The parameters are : 

1 ) inversion reduction factor gamma 

2) linear loss L 

3) initial transmission of absorber T[0] 

4) gain medium length I in cm 
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5) cavity length l_cav in cm 



6) beam radius in amplifier wO in cm 



>points := 
{seq([0.5+0.5*k/100,Pow(1.6,0.01,0.886,4.9,35,0.06,k)],k=1..100)}: 
>plot(points,x=0.5..1, style=point, symbol=circle, color=black, axes=boxed , 
title='Pulse power vs. R, [kW]'); 

Pulse power vs. R, [kW] 
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And, at last, the pulse durations are: 

>Width := proc(gam,L,T_0,l,l_cav,i) # definition of system's parameters 
>delta := 0.028: 

>sigma[g] := 7e-21:# in [cm 2 ] 

>sigma[s] := 3.5e-19:# Co:MALO 

> alpha := sigma[s]/(sigma[g]*gam): 

>h := 6.63e-34:# J*s 
>nu := evalf(3e8/1.354e-6):# frequency for 1.54 micrometers 
>S := Pi*w0 2 /2:# area for Gaussian beam 
>R := 0.5 + 0.5*i/100: 
>c := 3el0: 
>refractivity := 1.5:# active medium 
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>t_cav := 2*(l_cav-l)/c + 2*(l*refractivity)/c: 

>x_i := l/2*(ln(l/R) + L + ln(l/(T_0 2 )))/(l*sigma[g]): 

>x_tO := 
l/2*(ln(l/R) + L + delta*ln(l/(T_0 2 )))/(sigma[g]*l): 

>eql := 
-2*sigma[g]*x*l + ln(l/R) + L + (x/x_i) a *ln(l/(T_0 2 )) + delta*ln(l/(T_0 2 )) 
- delta*ln(l/(T_0 2 ))*(x/x_i) Q = 0:# e5 

>eq2 := 
ln(x)*delta*ln(l/(T_0 2 ))*alpha - ln(x_i)*delta*ln(l/(T_0 2 ))*alpha - 2*l*x* 
sigma[g]*alpha + 2*l*x_i*sigma[g]*alpha + ln(x)*L* alpha + ln(x)*ln(l/R)* 
alpha - ln(x_i)*L*alpha - ln(x_i)*ln(l/R)*alpha + (x/x_i) a *ln(l/(T_0 2 )) - 
delta* ln(l/(T_0 2 ))*(x/x_i) a - ln(l/(T_0 2 )) + delta*ln(l/(T_0 2 )) = 0:# e6 

>sol_f := fsolve(eq2, x, avoid={x=0}): 
>sol_t := fsolve(eql, x, avoid={x=0}): 
>sol_Width := 
(-l/2*ln(x_i/sol_f)*t_cav*alpha/(sigma[g]*l* (-x_i*alpha+sol_t*alpha+x_tO* 
ln(x_i/sol_t)*alpha + x_i - x_i*(sol_t/x_i) a - x_tO + x_tO* (sol_t/x_i) a )))* 
le9:# [ns] 

>end: 

>print('The parameters are:'); 
>print('l) reducing parameter gamma'); 
>print('2) linear loss L'); 

>print('3) initial transmission of absorber T[0]'); 
>print('4) gain medium length 1 in cm'); 
>print('5) cavity length l_cav in cm'); 

The parameters are : 

1) reducing parameter gamma 

2) linear loss L 

3) initial transmission of absorber T[0] 

4) gain medium length I in cm 

5) cavity length l_cav in cm 
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>points := 
{seq([0.5+0.5*k/100,Width(1.9,0.04,0.89,4.9,35,k)],k=1..100)}: 
>points_exp := {[0.793,70], [0.88,70], [0.916,75]}: 

>plot (points, x=0. 5.. l,style=point,axes=boxed,symbol=circle,title='Pulse width 
vs. R, [ns]'): 

>plot(points_exp,x=0.5..1,style=point,axes=boxed,color=red): 
>display({%,%%},view=50..100); 

Pulse width vs. R, [ns] 
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The worse agreement with the experimental data for the pulse durations is 
caused by the deviation of the pulse shape from the Gaussian profile, which was 
used for the analytical estimations. More precise consideration is presented in 

i- 

The main advantage of the analytical model in question is the potential of the 
Q-switched laser optimization without any cumbersome numerical simulations. 
Let's slightly transform the above obtained expressions: 



>e7 := subs( x=x_t,expand( e5/(2*sigma[g]*x_i*l) ) ); 
>e8 := subs( x=x_f,expand( rhs(sol_2)*l_cav*gamma/l 
>sol_6; 
>sol 3; 



= 0; 
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a;_i <jgX_il a g x_il o g x_il a g x_il 

2 <Jg x _il 



= 



-(51n(-^)ln(a;_/) £ln(— 7 )ln(*_0 -ln(4)ln(a; /) 
-a;_/ + ;r_» + - ^ i Jo , 2 R _ 



I <Tg 2 I (Jg I <Tg 

1 1 1 , , 1 s ( %_f ^C 

1 ln(-)ln(i_i) -Lln(i_/) 1 £ln ( T _ t ) 2 l T 2M x_ t J 

2 Z (7g Z <7g 2 I (Jg I <T g £ 

, ln(-^)d(^=4)< Iln(^),5 1 ln(-^) 
1 V T 2 x_z ; 2 V 2 1 V 2/ 



= 



'^) + L + < 



^(Ij + L + Jln^) 
a; tO 



x i 



2 CTgZ 



2 ^~] 



># Hence 

>e9 := x_t/x_i = (ln(l/R) + L + delta*ln(l/(T[0] 2 )))/ (2*l*x_i*sigma[g]) 
+ (x_t/x_i)«*ln(l/T[0] 2 )*(l-delta)/ (2*l*x_i*sigma[g]); 

>elO := x_f - x_i = m(x_f/x_i)*(ln(l/R) + L + delta*ln(l/(T[0] 2 )))/ 
(2*l*sigma[g]) - ln(l/T[0] 2 )*(l-delta)*(l-(x_f/x_i) c )/ (2*l*sigma[g]*alpha); 
> simplify ( sol_6 - sol_3 ); 

>ell := subs( 
ln(l/T[0] 2 ) = (x_i-x_tO)*2*l*sigma[g]/(l-delta),subs(ln(l/R) + L + delta* 
ln(l/(T[0] 2 )) = x_t0*2*l*sigma[g],e9)); 

>el2 := subs( 
ln(l/T[0] 2 ) = (x_i-x_t0)*2*l*sigma[g]/ (l-delta),subs(ln(l/R) + L + delta* 
ln(l/(T[0] 2 )) = x_t0*2*l*sigma[g],el0)); 



cj '.= — — — = — -+- — 

x_i 2 a g x_il (jgX_il 
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elO := x _f — x _i 

i ln( rr } (ln( fi } + L + Sln ^ i ln ^ (* - *) a - ( 



■''- 7 -' 1 ' ^ ' ' —-^)) ^n^Hl - S)(l - Mfr 

Jo i ^0 



I <7g 2 I (Jg a 



ln(-L)(-l + *) 

x W -x i " 



2 Z 



fa 



a; £ a; w x z — — 

ell := -^r = -^— H = : 



a; i a; i 



(»_» - x_tO) (1 - (^4) f ) 
el2 := x_f - x_i = ln(^=4) z_tt? — 



xi a 



Now, let's plot some typical dependence for Q-switching. 

>eq := subs({x_tO=x,x_i=l},subs(x_t/x_i=y,expand(ell)));#herex=x_tO/x_i, 

y=x_t/x_i 



eq := y = x + y^ — y^ x 



>xt := proc(zeta,i) 
>sol := array(1..20): 

>x := (i-l)/20: 

>sol[i] := fsolve(y = x+y^-y^*x,y,avoid={y=l},0.. infinity): 

>end: 

>points := {seq([k/20,xt(1.5,k)],k=2..19)}: 

>pl := plot(points,x=0..1,style=point,axes=boxed): 

>points := {seq([k/20,xt(3,k)],k=2..19)}: 

>p2 := plot(points,x=0..1,style=point,axes=boxed): 

>points := {seq([k/20,xt(6,k)],k=2..19)}: 
>p3 := plot(points,x=0..1,style=point,axes=boxed): 

>display({pl,p2,p3},view=0..1,title='x_t/x_ivs x_tO/x_ifor different zeta'); 
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x_t/x_i vs x_tO/x_i for different zeta 
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From the definition of x_t0 this dependence for growing £ tends to the linear 
one, which correspond to the maximal efficiency of the population inversion 
utilization. For the fixed absorption and emission cross-sections, the C parameter 
increases as a result of the laser beam focusing in the saturable absorber. Then 
( (^ oo), we have: 

>el3 := lhs(ell) = op(l,rhs(ell)); 

>e!4 := expand(lhs(e!2)/x_i) = expand(op(l,rhs(e!2))/x_i); 



el3 



x t x tO 



ln(M0z W 

1, X -J i x % z 

el4 := r — 1 = = — : 



and the output energy optimization can be realized by this simple way: 

># Energy optimization 

>el5 := x_f/x_i-l = ln(x_f/x_i)*rhs(sol_6)/x_i;# from el4 and expres- 
sion for x_t0 

>el6 := lhs(%) = ln(x_f/x_i)*(a + b + delta*c);# a = ln(l/R)/ 
(2*sigma[g]*x_i*l), b = L/(2*sigma[g]*x_i*l), c = ln(l/T[0] 2 )/ 
(2*sigma[g]*x_i*l) are the relative shares of the output, linear and saturable 
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loss in the net-loss 

>expand(solve(%,a)): 

>sol_7 := (y/log(y) - l/log(y) - (l-delta)*b - delta) /(l-delta);# solution 
for a, y = x_f/x_i, we used a+b+c=l 

>epsilon := -a*log(y);# dimensionless energy, 
epsilon = En*2*sigma[g]*gamma/(h*nu*A)/(2*sigma[g]*x_i*l) 
>simplify( subs(a=sol_7, epsilon) ); 
>diff(%,y) = 0;# maximum of energy 
> print ('An optimal ratio of final and initial inversions:'); 
>y_opt := solve(%,y);# optimal y 
>print('An optimal mirror:'); 

>a_opt := subs(y=y_opt,sol_7);# optimal a 
>print('A maximal dimensionless energy:'); 
>epsilon_opt := simplify ( subs({a=a_opt,y=y_opt}, epsilon) );# optimal ep- 
silon 

xf ^^(H^ + L + SM^)) 

elo := — — I = = 

X _l 2 <7glX_l 



el6 := ^=£ - 1 = In(— ) (a + b + Sc) 

x i x i 



so i 7:= Hy)_Hy) 



(l-S)b-S 



1-6 

e := — aln(y) 

y — 1 — bln(y) + bhi(y) S — Shi(y) 
-1 + 6 



b b6 6 

1- - + 

V V V 
-1 + 6 



= 



An optimal ratio of final and initial inversions 

y_opt := b — b 6 + 6 
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An optimal mirror : 



b-b5 + S 



ln(b-bS + 5) ln(b-bS + 5) 
a_opt := 



(l-S)b-S 



1-5 



A maximal dimensionless energy : 



epsilon _opt := 
b-bS + 6-l-b ln(6 - b 5 + 5) + b ln(6 -bS + 6)6-6 ln(6 - 6 (5 + 6) 



>pl := plot( 
{[b,subs(delta=0,a_opt),b=0..1],[b,l-subs(delta=0,a_opt)-b,b=0..1] 
},color=black,axes=boxed,title='optimal a and c versus b'): 

>p2 := 
plot( 

{[b,subs(delta=0.05,a_opt),b=0..1],[b,l-subs(delta=0.05,a_opt)-b, b=0. .1]}, 
color=red,axes=boxed,title='optimal a and c versus b'): 

>p3 := 
plot( {[b,subs(delta=0.1,a_opt),b=0..1],[b,l-subs(delta=0.1,a_opt)-b, b=0..1 ]}, 
color=blue,axes=boxed, title= 'optimal a and c versus b'): 

>display(pl,p2,p3,view=0..1); 

>pl := plot([b,subs(delta=0,epsilon_opt),b=0..1], 
axes=boxed,color=black, title='optimal epsilon versus b'): 

>p2 := plot([b,subs(delta=0.05,epsilon_opt),b=0..1], 
axes=boxed,color=red, title='optimal epsilon versus b'): 

>p3 := plot([b,subs(delta=0.1,epsilon_opt),b=0..1], 
axes=boxed,color=blue, title= 'optimal epsilon versus b'): 

>display(pl,p2,p3); 
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optimal a and c versus b 
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optimal epsilon versus b 
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From the first figure the optimization can be performed by the 
simple graphical way. We have to define the appropriate for our 
scheme laser's net-loss and to determine (measure or calculate) 
the intracavity linear loss. This gives the value of b, which is 
the relation of the linear loss to the net-loss. The upper group 
of curves gives the value of c, the lower curves give a (for the 
different contribution of the excited-state absorption, i. e. S). 
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2.3 Two-color pulsing 

Now let consider a more complicated situation, which corresponds to the two- 
color Q-switching due to presence of the stimulated Raman scattering in the 
active medium (for example, Yb 3+ :KGd(W04)2, see Q). In this case the an- 
alytical modelling is not possible, but we can use the numerical capacities of 
Maple. 

Let's the gain medium length is l g and the Raman gain coefficient is g. We 
assume the exact Raman resonance and neglect the phase and group-velocity 
effects. Then the evolutional equations for the laser and scattered intensities I p 
and 7 S , respectively, are: 

>restart: 
>with(DEtools): 
>with(plots): 

>eql := diff(In[p](z),z) = -g*In[s](z)*In[p](z):# evolution of the laser inten- 
sity 

>eq2 := diff(In[s](z),z) = g*In[p](z)*In[s](z):# evolution of the Stokes com- 
ponent intensity 

>sys := {eql,eq2}; 

>IC := {In[p](0)=In[p,0],In[s](0)=In[s,0]};# initial conditions 

>sol := dsolve(sys union IC,{In[s](z),In[p](z)}): 



d d 

sys ■= {—In s (z) = gln s (z)ln p (z), — In p (z) = -g In s (z) InJz)} 
az oz 



IC := {Jn 8 (0) = ln.,0, ln p (0) = In p . Q } 



The integration produces: 

>subs(sol,In[s](z)):# z is the distance in the crystal 
>sol_l := In[s](z) = simplify(%); 
>subs(sol,In[p] (z) ) : 

>sol_2 := In[p](z) = simplify (%); 



sol_l := In s (z) = 



sol_2 := Irip(z) 



In p , o + e( z 9 ( In p. °+ /n *. ")) In S) 
(In Pt0 + In a ,o) In P ,o 



In P:0 + e( z 3( In p.o+ins,o)) j Us0 
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There are the amplification of the scattered field and the depletion of the 
laser field, which plays a role of the pump for the Stokes component. 
Now let's transit from the intensities to the photon densities: 

>#phi is the photon densities, the laser component at 1033 nm, the Stokes 
component at 1139 nm 

>sol_3 := simplify ( subs( 
{z=2*l[g],In[p,0]=c*h*nu[p]*phi[p,0]/2, 
In[s,0]=c*h*nu[s]*phi[s,0]/2,In[s](z)=c*h*nu[s]*phi[s](z)/2}, 
2*sol_l))/(c*h*nu[s]); 

>sol_4 := simplify ( subs( 
{z=2*l[g],In[p,0]=c*h*nu[p]*phi[p,0]/2,In[s,0]=c*h*nu[s]*phi[s,0]/2, 
In[p](z)=c*h*nu[p]*phi[p](z)/2},2*sol_2) )/(c*h*nu[p]); 

,_ e('»g c M%i+^.^..°))^ (%i + I/8 ^ ) 
%1 := v p (f>p,o 

, , , / n {V p (j) p0 + V a (j) s0 )(j) p0 

sol_4 :— 4>p(z) = 



Vp 4> p , o + e(*» 9 c h ("" *p. 0+ " s ^.°^ v s (j) St0 



Hence the photon densities evolution due to the stimulated Raman scattering 
obeys: 

>eq3 := diff(phi[p](t),t)[raman] = ( subs( 
{phi[p,0]=phi[p](t),phi[s,0]=phi[s](t)},rhs(sol_4) ) - phi[p](t))/t[cav];# t[cav] is 
the laser cavity period 

>eq4 := diff(phi[s](t),t)[raman] = ( subs( 
{phi[p,0]=phi[p](t),phi[s,0]=phi[s](t)},rhs(sol_3) )- phi[s](t) )/t[cav]; 

(Vp<t>p(t) +V s (j) s {t))(j)p(t) 



3 __ Vp<t) p {t) + e ( ( 99 c ' l ( !/ p^(*)+^^(*)))i/ s (/) s (t) 

eq3 :— {— (pp(t)j 



Mt) 



OZ t cav 

e ( l a ach(%i+u.4>.W)) < i, s (t)(%l + i /s <l> a (t)) 

%l + e (' g gcM%i+^0 3 (t)))^^ft) 

tcav 

%l:=v p 4> p {t) 

The inverse and ground-state populations for the gain medium and the ab- 
sorber (Cr 4+ :YAG) obey (see previous subsection): 
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>eq5 := diff(n(t),t) = 
-gamma*sigma[g] *c*phi[p] (t) *n(t) ; 

>eq6 := diff(n[0](t),t) = -rho*sigma[a]*c*phi[p](t)*n[0](t);# rho=S[g]/S[a] 
is the ratio of the beam area in the gain medium to that in the absorber 



d_ 

at 



eq5 := ^n(i) = —ycr g c<f) p (t)n(t) 



d 
eq6 := ^7«o(i) = - P o c <f> p (t) n (t) 



Contribution of gain, saturable, output and linear intracavity loss to the 
field's evolution results in: 

>eq7 := diff(phi[p](t),t)[gain] = 
2* (sigma[g] *n(t) *l[g] - sigma[a] *n[0] (t) *l[a]) *phi[p] (t) /t [cav] ; 

>eq8 := diff(phi[p](t),t) [linear] = (-ln(l/R[p]) - L[p])*phi[p](t)/t[cav]; 

>eq9 := diff(phi[s](t),t) [linear] = 
(-ln(l/R[s]) - L[s] - 2*kappa*sigma[a]*n[0](t)*l[a])*phi[s](t)/t[cav]; 



« / d ,„ Q g n(t) l g - a a n (t) l a ) <j> p (t) 

6q := ^ VP^'ha-in = 2 



d (-H^-)-L P )Mt) 

eq8 := (— (j> p (t)) linear = 



\ r\ . f u \ " i i it/teat , 

at t 



cav 



d (-ln(—)-L s -2K<j a n (t)l a )(j) s (t) 



==(«: *.(*)) 



linear 



Oi rsv / / vinear , 

Ol I cav 

where L p and L s are the linear loss for the laser and Stokes components, re- 
spectively, R p and R s are the output mirror reflectivity at the laser and Stokes 
wavelengths, respectively, k is the reduction factor taking into account the de- 
crease of the loss cross-section at the Stokes wavelength relatively to the lasing 
one. 

Hence, the field densities evaluate by virtue of: 

>eqlO := diff(phi[p](t),t) = rhs(eq3) + rhs(eq7) + rhs(eq8); 
>eqll := diff(phi[s](t),t) = rhs(eq4) + rhs(eq9); 
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(v p cj> p {t) + v s cf> s (t))(j) p (t) 



eqlO:=^<t> p (t) 



., , m __ v p <j) p (t) +e ( - l s9ch(u P 4, p (t)+u s 4> s (t))) 1 , s( j )s ^ 



*""■' t 



cav 



2(<j g n(t)l g - cr a n (t) l a ) <p p (t) R p 



v cav v cav 

e (l g gch(%l+u s *.(*))) ^(t) (%1 + y s a ( t )) 



%l + e (I 9 flcfc(%l+ V . *.(*))) v ^ (j) 

egiJ := ^ </> s (£) = ^_ 



&(*) 



(-ln(— ) - L s - 2K(j a n n (t)l a )(j) s {t) 

1 £*s 

t cav 

%l:=v p (j) p {t) 

In the agreement with the results of the previous subsection we can transform 
the equation for the populations in the gain and absorption media: 

>diff(n[0](n),n) = subs( 
{n(t)=n,n[0](t)=n[0](n)},rhs(eq6)/rhs(eq5) ); 

>diff(n[0](n),n) = zeta*rho*n[0](n)/n;# zeta=sigma[a]/sigma[g]/gamma 
>dsolve({%},n[0](n)); 



d , . pa a n Q (n) 



dn 7 a g n 



d (pn {n) 

7-noW = 

on n 



n (n) = _(7in (c ' ,) 
>eql2 := n[0] = n[0,i]*(n/n[i]) c *' 5 ; 

eql2 := n = n .i (— ) 



n >(Cp) 



71 j 



>simplify( subs( phi[s](t)=0,subs( 
{n(t)=n[i],n[0](t)=n[0,i],phi[p](t)=0},expand( 
rhs(eqlO)/phi[p](t) ) ) ) ) = 0:# the generation start 

>numer(lhs(%)) = 0: 
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>eql3 := n[i] = subs( 
2*sigma[a]*n[0,i]*l[a]=ln(l/T[0] 2 ),solve(%,n[i]) ); 



1 
ln(- 

1 Vr 
eql3 := n 



2 a a L 



1 T R P 



'g u g 



Let' introduce new variables: 

>Phi[p](t) = 2*sigma[g]*l[g]*phi[p];# note, that Phi is not the dimensional 
intensity from the first subsection! 
>Phi[s](t) = 2*sigma[g]*l[g]*phi[s]; 
>tau = 2*sigma[g]*l[g]*n[i]*t/t[cav]; 
>Nu = nu[s]/nu[p]; 

>G = g*l[g]*c*h*nu[p]/(2*sigma[g]*l[g]); 
>y(tau) = n(t)/n[i]; 
>Xi = ln(l/(T[0] 2 ))+ln(l/R[p])+L[p]; 
kappa = 0.38; 



*p(*) = 
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k = .38 



Then from Eqs. (eq5, eqlO, eqll, eql2): 

>f[l] := diff(Phi[p](tau),tau) = Phi[p](tau)/Xi* 
(Xi*y(tau) - ln(l/T[0] 2 )*y(tau) a *" - (ln(l/R[p])+L[p]) + 
( (Phi[p] (tau) +N*Phi[s] (tau) ) / (Phi[p] (tau) + 
exp(G*(Phi[p](tau)+N*Phi[s](tau)))*N*Phi[s](tau)) - 1 ) ); 

>f[2] := diff(Phi[s](tau),tau) = 
Phi[s](tau)/Xi*(-(ln(l/R[s])+L[s]+kappa*ln(l/(T[0] 2 ))* 
y(tau) a * p ) + 

(Phi[p] (tau) +N*Phi[s] (tau) ) / (Phi[p] (tau) * 
exp(-G*(Phi[p](tau)+N*Phi[s](tau))) + N*Phi[s](tau) ) - 1 ); 

>f[3] := diff(y(tau),tau) = -gamma*(l[cav]/l[g])*Phi[p](tau)/Xi; 



Lu,^( a ,)_wl. r $ p (t)+JY*.(t) 



$ p (r) (sy(r) - ln(-^) y(r)^ - ln(-l) - L 



*p(t) + e( G (*»M+w *.M)) TV $ s (r) 



/ 2 :=|:$ s (r) 

1 , , ..!,„■ ! u. M (<«p) , <IV(v) + .V<l»4r) 



$ s (r) -ln( — )-L s -«ln(-^)y(r)( Q ^ + 



i?/ s v T 2; ' rv y ^(rJeC-^^W+^^W^+JV*.^) 



<9t La 



The next procedures solve the obtained system numerically to obtain the de- 
pendencies of the normalized photon densities at laser and Stokes wavelengths 
and the relative inversion vs. cavity roundtrip, n is the duration of the simula- 
tion in the cavity roundtrips): 

>ODE_plotl := proc(TO,Rp,Rs,Lp,Ls,n) 
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>gam := 1:# gamma parameter 

>rho := 1: 

>kappa := 0.38: 

>alpha := evalf( 5e-18/2.8e-20 ):#the ratio of the cross-sections 

>x := 29:# x = l[cav]/l[g] 

>N := evalf(1033/1139): 

>G := subs( 
{g=4.8e-9,c=3el0,h=6.62e-34,nu_p=3el0/1.033e-4, 
sigma_g=2.8e-20},g*c*h*nu_p/(2*sigma_g) ): 

>Xi := ln(l/(T0 2 ))+ln(l/Rp)+Lp: 

>sys := [ D(Phip)(tau) = 
Pmp(tau)*(Xi*y(tau)-ln(l/(T0 2 ))*y(tau) Q *'' 
-ln(l/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau) + 
exp(G*(Pmp(tau)+N*Pms(tau)))*N*Pms(tau))-l)/Xi, 

D(Phis)(tau) = 
PmsttauHdnCl/Rsj-Ls-^ppa^ntl/tTO^^yttau)"*'' 
+ (Phip(tau)+N*Phis(tau))/(Phip(tau)* 
exp(-G*(Pmp(tau)+N*Pms(tau)))+N*Phis(tau))-l)/Xi, 

D(y)(tau) = -gam*x*Phip(tau)/Xi]: 

>DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n, 
[[Phip(0) = le-7,Phis(0) = le-7,y(0) = l]],stepsize=0.1, 
scene=[tau,Phip(tau)],method=classical[abmoulton], 
axes=FRAME,linecolor=BLUE) : 

>end: 

>ODE_plot2 := proc(TO,Rp,Rs,Lp,Ls,n) 

>gam := 1: 

>rho := 1: 

>kappa := 0.38: 

>alpha := evalf( 5e-18/2.8e-20 ): 

>x := 29:# x = l[cav]/l[g] 

>N := evalf(1033/1139): 

>G := subs( 
{g=4.8e-9,c=3el0,h=6.62e-34,nu_p=3el0/1.033e-4,sigma_g=2.8e-20}, 
g* c *h*nu_p/(2*sigma_g)): 

>Xi := ln(l/(T0 2 ))+ln(l/Rp)+Lp: 

>sys := [ D(Phip)(tau) = 
Phip(tau)*(Xi*y(tau)-ln(l/(T0 2 ))*y(tau) Q *''- 
ln(l/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau) + 
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exp(G*(Phip(tau)+N*Phis(tau)))*N*Phis(tau))-l)/Xi, 

D(Phis)(tau) = 
Phis(tau)*(-ln(l/Rs)-Ls-kappa*ln(l/(T0 2 ))*y(tau) Q *''+ 
(Phip(tau)+N*Phis(tau))/(Phip(tau)* 
cxp(-G*(Phip(tau)+N*Phis(tau)))+N*Phis(tau))-l)/Xi, 

D(y)(tau) = -gam*x*Phip(tau)/Xi]: 

>DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n, 
[[Phip(0) = le-7,Phis{0) = le-7,y(0) = l]],stepsize=0.1, 
scene=[tau,Phis(tau)],method=classical[abmoulton], 
axes=FRAME,linecolor=RED) : 

>end: 

>0DE_plot3 := proc(TO,Rp,Rs,Lp,Ls,n) 

>gam := 1: 

>rho := 1: 

>kappa := 0.38: 

>alpha := evalf( 5e-18/2.8e-20 ): 

>x := 29:# x = l[cav]/l[g] 

>N := evalf(1033/1139): 

>G := subs( 
{g=4.8e-9,c=3el0,h=6.62e-34,nu_p=3el0/1.033e-4,sigma_g=2.8e-20}, 
g*c*h*nu_p/(2*sigma_g)): 

>Xi := ln(l/(T0 2 ))+ln(l/Rp)+Lp: 

sys := [ D(Phip)(tau) = 
Phip(tau)*(Xi*y(tau)-ln(l/(T0 2 ))*y(tau) Q *'' 
-ln(l/Rp)-Lp+(Phip(tau)+N*Phis(tau))/(Phip(tau) + 
exp(G*(Phip(tau)+N*Phis(tau)))*N*Phis(tau))-l)/Xi, 

D(Phis)(tau) = 
Phis(tau)*(-ln(l/Rs)-Ls-kappa*ln(l/(T0 2 ))*y(tau) Q *'' 
+ (Phip(tau)+N*Phis(tau))/(Phip(tau)* 
exp(-G*(Phip(tau)+N*Phis(tau)))+N*Phis(tau))-l)/Xi, 

D(y)(tau) = -gam*x*Phip(tau)/Xi]: 

>DEplot(sys,[Phip(tau),Phis(tau),y(tau)],tau=0..n, 
[[Phip(0) = le-7,Phis(0) = le-7,j/(0) = l]],stepsize=0.1, 
scene=[tau,y(tau)],method=classical[abmoulton], 
axes=FRAME,linecolor=BL ACK) : 

>end: 
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>display(ODE_plotl(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed, 
title='laser photon density'); 

>display(ODE_plot2(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed, 
title='Stokes photon density'); 

>display(ODE_plot3(0.65,0.75,0.96,5e-2,5e-2,300),axes=boxed, 
title= 'relative inversion'); 
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So, we can obtain a quite efficient two-color pulsing in the 
nanosecond time domain without any additional wavelength 
conversion. 



The nonstationary lasing in question produces the pulses with the durations, 
which are larger than the cavity period. The cavity length shortening, i. e. 
the use of the microchip lasers, decreases the pulse durations down to ten - 
hundred picoseconds. But there is a method allowing the fundamental pulse 
width reduction, viz. mode locking. 

3 Conception of mode locking 



The laser cavity is, in fact, interferometer, which supports the propagation of 
only defined light waves. Let consider a plane wave, which is reflected from a 
laser mirror. The initial wave is (cc is the complex conjugated term): 

>restart: 

with (plots): 

with(DEtools): 

AI := A0*exp(I*(omega*t+k*x))+cc; 



AI := A0e {I{ut+kx » + cc 
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Then the reflected wave (normal incidence and full reflectivity are supposed) 

is: 

>AR := AO*exp(I*(omega*t-k*x+Pi))+cc; 

AR:= A0e {I{uJt - kx+7r)) + cc 
where 7r is the phase shift due to reflection. An interference between incident 
and reflected waves results in 



>convert(AI+AR,trig): 
expand (%): 
factor(%); 



-2A0 sin(w t) sin(fc x) + 21 AO cos(w t) sin(fc x) + 2 cc 



This is the so-called standing wave: 

>animate(sin(x)*sin(t), x=0..2*Pi,t=0..2*Pi, axes=boxed, 
title='Standing wave', color=red); 



Standing wave 
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-0.5 
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Note, that a wave node lies on a surface (point x=0). The similar situation 
takes place in the laser resonator. But the laser resonator consists of two (or 
more) mirrors and the standing wave is formed due to reflection from the each 
mirror. So, the wave in the resonator is the standing wave with the nodes placed 
on the mirrors. Such waves are called as the longitudinal laser modes. The laser 
resonator can contain a lot of modes with the different frequencies (but its nodes 
have to lie on the mirrors!) and these modes can interfere. 

Let suppose, that the longitudinal modes are numbered by the index m. In 
fact, we have M harmonic oscillators with the phase and frequency differences 
dphi and domega, correspondingly. Let the amplitude of modes is AO. 

>mode := l/2*A0*exp(I*(phi0+m*dphi) + 

I*(omegaO+m*domega)*t)+cc;# amplitude of the mode numbered 
by index m 



mode ■— — AO pi 1 (4>°+ m d P hi )+ I M+™ domega) t) , cc 



Here 4>0 and w0 are the phase and frequency of the central mode, respectively. 
The interference between these modes produces the wave packet: 

>packet := sum(mode-cc,m=-(M-l)/2..(M-l)/2)+cc; 

# interference of longitudinal modes with constant phases 



1 AO e^ ) e^ *) e( / ^ 1 / 2M+1 / 2 ^ dphi ^ e ( It ( 1 / 2M + 1 / 2 ) domega) 

" ' O p(I dphi) pil domega t) 1 

1 AO e( 7< ^°) e ( 7w0t ) e ( / (- 1 / 2M + 1 / 2 ) d P hi ) e (It (-1/2 M+l/2) domega) 
O p(I dphi) P (I domega t) 1 



CC 



Now, we can extract the term describing the fast oscillation on the central 
("carrier") frequency wO from the previous expression. The obtained result is 
the packet's envelope (its slowly varying amplitude): 

>envelope := expand((packet-cc)/exp(I*(t*omegaO+phiO))); 
# slowly varying envelope of the wave packet 



envelope := — 



1 AO e^ 1 / 21 d P hiM ) p^-I 21 dphi) g(l/2i domegatM) e (l/2 J domegat) 



2 p(I dphi) p(I domegat) 1 

1 AO e^ 1 / 27 dphiM) g(l/2 / dphi) J-1/2I domegat M) g(l/2 7 domegat) 

2 p (I dphi) p(l domegat) 1 
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It is obviously, that this expression can be converted into following form: 

>envelope := 

l/2*A0*sinh(l/2*I*M*(dphi+t*domega))/ 
sinh(l/2*I*(dphi+t*domega)); 



envelope := — 
2 



-i AO sin(— M (dphi + domegat)) 



sin(— dphi -\ — domegat) 



The squared envelope's amplitude (i. e. a field intensity) is depicted in the next 
figure for the different M. 

>plot({subs({M=5,A0=l,dphi=0.1,domega=0.1},2*envelope 2 ), 
subs({M=20,A0=l,dphi=0.1,domega=0.1},2*envelope 2 )}, 
t=-80..80, axes=boxed, title='result of modes interference'); 
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result of modes interference 
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One can see, that the interference of modes results in the generation of short 
pulses. The interval between pulses is equal to 2 n/domega. The growth of 
M decreases the pulse duration 2 ir/(M*domega) and to increases the pulse 
intensity M 2 *A0 . The last is the consequence of the following relation: 



> 



(limit(sin(M/2*x)/sin(x/2),x=0)*A0) 2 ;# maximal field intensity 
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M 2 A0 2 
In this example, the phase difference between neighboring modes is constant. 

Such mode locking causes the generation of short and intense pulses. But in the 
reality, the laser modes are not locked, i. e. the modes are the oscillations with 
the independent and accidental phases. In this case: 

>M := 20:# 20 longitudinal modes 

A0 := 1:# constant dimensionless amplitude of mode 

domega := 0.1 :# constant frequency different 
phiO := 0:# phase of the central mode 
omegaO := 1:# dimensionless frequency of central mode 
mode := 0: 

for m from -(M-l)/2 to (M-l)/2 do: 
die := rand(6): 
dphi := die():# accidental phase difference between modes 

mode := mode+A0*cos(phi0+m*dphi+(omega0+m*domega)*t): 

od: 
plot(mode,t=-80..80, axes=boxed, title='result of modes interference'); 



result of modes interference 




-60 -40 -20 



Thus, the interference of the unlocked modes produces the irregular field beat- 
ings, i. e. the noise spikes with a duration ~ l/(M*domega). 
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What are the methods for the mode locking? Firstly, let consider the sim- 
plest model of the harmonic oscillation in the presence of the periodical force. 

>diff_equation := diff(diff(y(t),t),t)+omega 2 *y(t)=cos(delta*t+phi); 

# oscillations in the presence of periodical force 

# (delta and phi are the frequency and 
#phase of the force oscillation, respectively) 

dsolve({diff_equation, y(0)=l, D(y)(0)=0},y(t)): 
oscilll := combine (%); 



diff _equation := (tt-jt y(t)) + u) y(t) = cos(<5 1 + ■ 



oscilll := y(t) = (us cos(oj t — <fi) + oj cos(w t + <f) — 

2 cos(w t) a; 3 + 2 u cos(a; t) 5 2 — S cos(uj t — <j>) 
+ S cos(uj t + 4>)-2 cos(<5 t + 4>)uj) /(-2 uj 3 + 2 uj S 2 ) 

>animate(subs({omega=l,delta=0.1 },subs(oscilll,y(t))), 
t=0..100,phi=0..2*Pi,color=red,style=point, axes=boxed); 




100 



We can see, that the external force causes the oscillations with the additional 
frequencies: uj+ S and uj- 8. If 5 is equal to the intermode interval, the addi- 
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tional oscillation of mode plays a role of the resonance external force for the 
neighboring modes. Let consider such resonant oscillation in the presence of the 
resonant external force: 

>diff_equation := diff(diff(y(t),t),t)+omega 2 *y(t) = 
cos(omega*t+phi);# resonant external force 

dsolve({diff_equation,y(0)=l,D(y)(0)=0},y(t)): 
oscill2 := combine(%); 



d 2 
diff _ equation :— (—-^y(t))+uj y(t) — cos(uit + i 



oscill2 := y(t) = 
1 — cos(uj t — <j>) + cos(w t + <fi) + 4 cos(w t) lo 2 + 2 1 s'm(oj t + <f>) ut 
4 ~ 2 

The term, which is proportional to t ("secular" term), equalizes the phase of the 
oscillations to the phase of the external force. It is the simplest model of a so- 
called active mode locking. Here the role of the external force can be played by 
the external amplitude or phase modulator. Main condition for this modulator 
is the equality of the modulation frequency to the intermode interval that causes 
the resonant interaction between modes and, as consequence, the mode locking 
(part 4)- The different mechanism, a passive mode locking, is produced by 
the nonlinear interaction of modes with an optical medium. Such nonlinearity 
can be caused by saturable absorption, self-focusing etc. (see further parts of 
article). Now we shall consider the simplest model of the passive mode locking. 
Let suppose, that there are two modes, which oscillate with different phases and 
frequencies in the cubic nonlinear medium: 

>omega := 1: 
delta := 0.1: 

sys := [diff(diff(y(t),t),t)+omega 2 *y(t)=-(y(t) 3 +z(t) 3 ), 
diff(diff(z(t),t),t) + (omega+delta) 2 *z(t)=-(y(t) 3 +z(t) 3 )]; 
#two oscillating modes with nonlinear coupling 

DEplot3d(sys,[y(t),z(t)],t=0..100,[[y(0)=-l,z(0)=l,D(y)(0)=0, 

D (z) (0) =0]] ,stepsize= 1 ,scene= [y (t) ,z(t) ,t] , 

axes=boxed,linecolor=BLACK,orientation=[-60,70], 

title='mode locking'); 



sys := 

[(|i y(*)) + y(*) = -y W 3 - Z W 3 - (^ Z W) + 1.21 <t) = - y (t) 3 - z (t) 3 
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mode locking 




We can see, that the initial oscillations with the different phases are locked due 
to nonlinear interaction that produces the synchronous oscillations. 



As conclusion, we note that the mode locking is resulted from 
the interference between standing waves with constant and 
locked phases. Such interference forms a train of the ultrashort 
pulses. The mechanisms of the mode locking are the external 
modulation with frequency, which is equal to intermode fre- 
quency interval, or the nonlinear response of the optical medium 



Later on we shall consider both methods. But firstly we have to obtain the more 
realistic equations describing the ultrashort pulse generation. 



4 Basic model 



The models describing the laser field evolution are based usually on the so- 
called semi-classical approximation. In the framework of this approximation 
the field obeys the classical Maxwell equation and the medium evolution has 
the quantum-mechanical character. Here we shall consider the wave equation 
without concretization of the medium evolution. 
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The Maxwell equation for the light wave can be written as: 

>restart: 

with (PDEtools ,dchange) : 

maxwell _eq := 
diff(E(z,t),z,z)-dirT(E(z,t),t,t)/c 2 =4*Pi*diff(P(t),t,t)/c 2 ; 



maxwell eq := (^ E(z, t)) - m - \ =4 SM-^ll 

oz z c c z 



where E(z,t) is the field strength, P(%) is the medium polarization, z is the 
longitudinal coordinate, t is the time, c is the light velocity. The change of the 
variables z > z*, t - z/c -> t* produces 

>macro(zs='z*',ts='t*'): 

tr := {z = zs, t = ts + zs/c}; 
maxwell_m := dchange(tr,maxwell_eq,[zs,ts], simplify); 



z* 

tr := {t = t* H , z = z*} 

c 



>e// 



(Jg E(z*, f *, c))c-2 (g^ E(s*, f *, c)) 

maxwell m := -ozz 

c 

W^P^*, «*, c)) 



We does not take into account the effects connected with a wave propagation in 
thin medium layer, that allows to eliminate the second-order derivation on z*. 
Then 

>int(op(2,expand(lhs(maxwell_m))),ts)-int(rhs(maxwell_m),ts): 
# integration of both sides of wave equation 
numer(%): 

wave_l := expand(%/(-2)); 

d d 

wave 1 := (— — E(z*, t*, c)) c + 2 n (— — P(z*, £*, c)) 

az* c?£* 

So, we reduced the order of wave equation. The inverse transformation of the 
coordinates leads to the so-called shortened wave equation: 
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>tr := {zs = z, ts = t - z/c}; 

wave_2 := dchange(tr,wave_l,[z,t], simplify); 



r z i 

tr := {z* = z, t* = t } 

c 



wave_2 := (— E(z, t, c)) c + (— E(z, t, c)) + 2^ (— P(z, t, c)) 



Next step is the transition to the slowly-varying amplitude approximation. We 
shall consider field envelope p(z,t) and polarization P(t), which are filled by the 
fast oscillation with frequency ui (k is the wave number, N is the concentration 
of the atoms, d is the medium length): 

>Theta=omega*t-k*z;# phase 

E := rho(z,t)*cos(Theta);# field 
P := N*d*(a(t)*cos(Theta)-b(t)*sin(Theta));# polarization (a and b are the 
quadrature components) 



G = uot — k z 



E := p(z, f)cos(e) 



P:= Nd (a(i) cos(9) - b(t) sin(0)) 

Then the wave equation can be transformed as: 

>Theta:=omega*t-k*z: 
diff(E,z)+diff(E,t)/c+2*Pi*diff(P,t)/c: 
combine(%,trig): 

collect(%,cos(Theta)): 
wave_3 := collect(%,sin(Theta)): 
eq_field := 

op(l,expand(coeff(wave_3,cos(Theta))))+op(3,expand( 

coeff(wave_3,cos(Theta)))) = 
-op(2,expand(coeff(wave_3,cos(Theta))))- 
op(4,expand(coeff(wave_3,cos(Theta))));# field equation 

disp_conditions := Int(coeff(wave_3,sin(Theta)),t)=0;# dispersion con- 
dition 
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e q field := (#. ,(,, t)) + M^l = _ 2 -JlMam + 2^db( t )o, 

az c c c 



/ 



disp _ conditions := 
p(z, t)kc-2nNd{£ I b{t)) - p(z, t) uj - 2 tt N d a,(t) lu 



The obtained result is the system of the shortened wave equa- 
tion and the dispersion condition. The right-hand side of the 
field equation (material part) will be different for the different 
applications (see below) 



5 Active mode locking: harmonic oscillator 

5.1 Amplitude modulation 

We start our consideration with a relatively simple technique named as active 
mode locking due to amplitude modulation. The modulator, which is governed 
by external signal and changes the intracavity loss periodically, plays the role of 
the external "force" (see part 2). Let consider the situation, when the modula- 
tion period is equal and the pulse duration is much less than the cavity period. 
If the modulation curve is close to cosine, then the master equation for the 
ultrashort pulse evolution can be written as []5|: 

£p(z, t)=g p{z, t) + t f 2 j^p{z, t) -Alt 2 p(z, t), 

where g is the net-gain in the laser accounting for the saturated gain a and 
linear loss / (including output loss), tf is the inverse bandwidth of the spectral 
filter, M is the characteristic of the modulation strength taking into account 
curvature of the modulation curve at the point of maximal loss. 

Let try to solve this equation in the case of the steady-state propagation of 
ultrashort pulse (when J| p(z, i)=0) and in the absence of detuning of mod- 
ulation period from cavity period. If the time is normalized to tf , then the 
obtained equation is a well-known equation of harmonic oscillator: 

>restart: 

with (plots): 
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with('linalg'): 
ode[l] := diff(rho(t),'$'(t,2)) + g*rho(t) - M*t 2 *rho(t); 



n2 



>sol := dsolve(ode[l]=0,rho(t)); 



sol := p(t) 



CI WhittakerW( 1 - 


9 1 

v/M' 4' 




Mt 2 ) 


Vt 

_C2WhittakerM( 


1 9 

4 VM' 


1 

4' 


,VMt 2 ) 



Vt 

The next step is suggested by the asymptotic behavior of the prospective 
solution: lim^oo p(t) = 0. But previously it is convenient to transit from 
Whittaker functions to hypergeometric and Kummer functions: 

WhittakerM( fi, f,z)= e(~> x^ + ^ hypergeom( i + v — (A, 1+2 v, x), 
WhittakerW( /x, v,z) = e^"3") x(2 +l/ ) KummerU( | + f — //, 1+2 v, x). 

As result, we have ( p,= j= , v= j, x= \/Mt 2 ): 



p(t) = Cie(-**M i VVI7 hypergeom( f - -$_, |, VMi 



4 4v / M' 2' 



C 2 e^^)vVMKummerU( | - ^4— , |, VMt 2 ) = 
C x e(-^i/7M h yP ergeom(( ± - ^_)+ §, |, VMi 2 ) 
C 2e <-^)VVS hy P ergeom( | - -f=, i, VMt 2 ). 



Now, the asymptotic condition lim^oo p(t) = 0, which is similar to condition 
for quantum states in harmonic oscillator, results in (see, for example, 0) 



pi(t) = C HermiteH 2 „(t vW) e^^ 31 ) for n= - {\ - j3=) && n is 

integer, 



p 2 (t)= C HermiteH 2n+l (t \HM) e^Mr 11 ) for n= - ( § - ^4=) && n is 

integer, 
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where HermiteH2n and HermiteH 2n+i are Hermite polynomials. Value of 
constant C can be obtained from the energy balance condition, which results 
from the equation of gain saturation: 



~ l+X JZo P(*) 2 dt ■ 

Here % is the inverse flux of the gain saturation energy, ao is the gain for small 
signal defined by gain medium properties and pump intensity (note, that g = 
a — I). Now let investigate the parameters of the steady-state solution of the 
master equation. With that end in view, we have to search the generation field 
energy: 



>C 2 *int( HermiteH(0,t*surd(M, 4)) 2 * exp( -t 2 *sqrt(M) 
t=-infinity.. infinity ); # energy of the first solution 
("ground state"), where alpha=l+sqrt(M) (see above) 



C 



2 



Af(V4) 



The use of normalization of intensity p 2 to (x*/) an d energy balance con- 
dition (see above) results in 

>eq := l+sqrt(M) - alpha[0]/(l+C 2 *sqrt(Pi)/surd(M, 4)) = 0: 
sol := solve (eq, C 2 );# pulse peak intensity 



surd(M, 4) (-1 - VM + ao) 
sol := 



Hence the ultrashort pulse obtained as result of active mode locking can be 
represented as 

>pulse := sol*HermiteH(0,t*surd(M, 4)) 2 * exp( -t 2 *sqrt(M) ); 

# this is first solution of master equation, 

# which represents a Gaussian pulse profile 

plot( subs({l=0.1,M=0.05,alpha[0]=1.2},pulse), 
t=-6..6, axes=BOXED, title='field intensity vs time' ); 

surd(M, 4) (-/ - VM + a ) e^^^ 
pulse := 



/7r(/ + VM) 
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field intensity vs time 




Now we calculate the pulse duration measured on half-level of maximal pulse 
intensity: 

>eq := exp( -sqrt(M)*t 2 ) = 1/2: 
sol := solve(eq, t): 

pulse_width := simplify( sol[l] - sol[2] ); 



pulse _width := 2 



y/Af (3/ 2 ) ln(2) 
M 



>plot( pulse_width,M=0.3..0.01, axes=BOXED, 

title='pulse width vs M'); 
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pulse width vs M 




0.3 



One can see, that the increase of modulation parameter decreases the pulse 
width, but this decrease is slow (~1/ vVM)- Next step is the taking into 
account the detuning of cavity and modulation periods. The corresponding 
normalized parameter S can be introduced in the following form (see corre- 
sponding Doppler transformation: t->t - z S and, as consequence, ^ ->- 8 ^ 
for steady-state pulse, here S is, in fact, inverse relative velocity with dimension 
[time/cavity transit], i. e. time delay on the cavity round-trip): 



>ode[2] := diff(rho(t),T(t,2)) 
g*rho(t) -M*t 2 *rho(t); 



delta*diff(rho(t),t) 



o2 pi 

ode 2 :=(—p(t)) + 6(-p(t))+gp(t)-Mt 2 p(t) 



>sol := dsolve(ode[2]=0,rho(t)); 



CI WhittakerW(— 4g ** , -, >/M r * 2 )e {_1/2 " ) 



sol := p(t) 



C2 WhittakerM( 



•16 VM ' 4' 

Vt 

1 4g-6 2 1 
16 JM ' 4 



v^ 



,VMt 2 )e ( - 1/25t) 
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The comparison with above obtained result leads to 



pi (i) = C-HermiteH 2 „(iV / v / M)e ( 2 ) 



for 



1 ^-43, 
4 16-y/M 



&& n is integer, 



, / t (5+VMt) , 

p 2 (t) = C*-HcrmiteH2„ +i (f\/\/M)e ( § ) 



for 



,3 , 5 2 -4g , 

n = -(7 ^ 7=) 

4 16VM 



&& n is integer, 



and we can repeat our previous analysis 

>C 2 *int( HermiteH(0,t*surd(M, 4)) 2 * 

exp( -t*(delta+sqrt(M)*t) ),t=-infinity.. infinity 

# energy of the first solution ("ground state"), 

# where alpha=l+sqrt(M)+delta 2/4 



^,2 (1/4-L) /- 



M(V4) 

>eq := l+sqrt(M)+delta 2 / 4 - alpha[0]/(l+%) = 0: 
# energy balance condition 

sol_0 := solve(eq, C 2 );# pulse peak intensity 



sol 



M^ 1 / 4 ) (4 1 + 4 VM + (5 2 - 4 ao) 



,(1/4 
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V^(4/+4VM + <5 2 ) 
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Note, that there is the maximal \S\ permitting the "ground state" ultrashort 
pulse generation: 

>-4*l-4*sqrt(M)-delta 2 +4*alpha[0] > 0; 

solve(-4*l-4*sqrt(M)-delta 2 +4*alpha[0]=0, delta); 



0< -41-4VM -S 2 +4 



»o 



2 V -I - v M + a , -2\/-l- VM + a 



We see, that the upper permitted level of detuning parameter is increased due 
to pump growth (rise of ao) and is decreased by M growth. 

It is of interest, that the growth of |<5| does not leads to generation of the 
"excite state" pulses because of the corresponding limitation for these solutions 
is more strict: 

>En := C 2 *int( expand( HermiteH(l,t*surd(M, 4) )) 2 * 
exp(-t*(delta+sqrt(M)*t) ), t=-inflnity.. infinity ): 
# energy of the second solution ("excite state") 

eql := = - (3/4+(delta 2 -4*(alpha-l))/(16*sqrt(M))): 
solve (%, alpha): 

eq := % - alpha[0]/(l+En) = 0:# energy balance condition 
solve(eq, C 2 );# "excite state" pulse peak intensity 

-12*sqrt(M)-4*l-delta 2 +4*alpha[0] > 0;# condition of generation 



M< 5 / 4 ) (12 VM + 4 / + 5 2 - 4 a ) 



e^ 1 / 4 7^) SU rd(M, 4) 2 0F (14 VM S 2 + 24 M + 5 A + 4 1 5 2 + 8 1 VM) 



< -12VM-4/-(5 2 +4a 



The dependence of pulse width on S has following form: 

>eq := exp( -t*(delta+sqrt(M)*t) ) = 1/2: 

sol := solve(eq, t): 

pulse_width := simplify ( sol[l] - sol[2] ); 

plot( subs(M=0.05,pulse_width),delta=-5..5, axes=BOXED, 
title= 'pulse width vs delta' ); 
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pulse _width 



^M{S 2 +4\/Mln(2)) 



M 



pulse width vs delta 




We can see almost linear and symmetric rise of pulse width 
due to detuning increase. This is an important characteristic 
of active mode-locked lasers. But in practice (see below), the 
dependence of the pulse parameters on detuning has more com- 
plicated character 



The dependence of the pulse maximum location on detuning parameter can be 
obtained from the solution for pulse envelope's maximum: -£■ e { 2 > — Q. 



Hence, the pulse maximum location is: 



>diff(exp(-t*(delta+sqrt(M)*t)/2), t) = 0: 
solve(%, t); 



1_5 
2 



>plot( subs(M=0.05,%),delta=-5..5, axes=BOXED, 
title= 'pulse shift vs delta' ); 
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pulse shift vs delta 




The increase (decrease) of the pulse round-trip frequency ( 5<0 and <5>0, respec- 
tively) increases positive (negative) time shift of the pulse maximum relatively 
modulation curve extremum. 



5.2 Phase modulation 



The external modulation can change not only field amplitude, but its phase. 
In fact, this regime (phase active modulation) causes the Doppler frequency 
shift of all field components with the exception of those, which are located in 
the vicinity of extremum of modulation curve. Hence, there exists the steady- 
state generation only for field located in vicinity of points, where the phase is 
stationary. Steady-state regime is described by equation 



>ode[3] := diff(rho(t),T(t,2)) 
I*M*t 2 *rho(t); 



g*rho(t) + I*phi*rho(t) 



d 2 

ode 3 := (-^ p{t)) + g p(t) + 1 p{t) + I M t 2 p(t) 



where <f> is the phase delay on the cavity round-trip. This equation looks like 
previous one, but has complex character. This suggests to search its partial 
solution in the form p(t) = Ce ((a+;l,)t!) : 
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>expand( subs(rho(t)=C*exp((a+I*b)*t 2 ), ode[3]) ): 
expand(%/C/exp(t 2 *a)/exp(I*t 2 *b)): 
eql := collect (coeff(%,I), t 2 ): 

eq2 := collect ( coeff(%%,I,0), t 2 ): 
eq3 := coeff(eql, t 2 ); 
eq4 := coeff(eql, t,0); 
eq5 := coeff(eq2, t 2 ); 
eq6 := coeff(eq2, t,0); 

soil := solve(eq6=0, a);# solution for a 

sol2 := solve(subs(a=soll,eq3=0), b);# solution for b 
sol3 := solve(subs(b=sol2,eq4)=0, phi);# solution for phi 
# but NB that there is eq5 defining, in fact, g: 
solve( subs({a=soll, b=sol2},eq5)=0, g);# solution for g 



eq3 :=8ab + M 

eq4 :=2b + 4> 
eq5 :=4a 2 -4b 2 



eq6 := 2 a 



soil := a 

2 y 



1 M 
sol2 := 

4 9 



1 M 
solS ----- 

2 9 



- V-2M, -- V-2M, - V2~VM, -- V2~VM 



So, we have p(t) = CeW / / and </> = — w — ,5= 

^ . The time-dependent parabolic phase of ultrashort pulse 
is called chirp and is caused by phase modulation 
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Pulse width for this pulse is: 

>exp(-sqrt(M/2)*t 2 ) = 1/2: 
sol := solve (%, t): 

pulse_width := simplify( sol[l] - sol[2] 



pulse _width 



^M^/ 2 ) V21n(2) 



M 



that is 2 i /ln(2) J jj . As one can see, that result is equal to one for amplitude 

modulation. The main difference is the appearance of the chirp. 

Let take into account the modulation detuning. In this case we may to 
suppose the modification of the steady-state solution: 

>ode[4] := diff(rho(t),T(t,2)) + g*rho(t) + delta*diff(rho(t),t) + 
I*phi*rho(t) + I*M*t 2 *rho(t); 

expand( subs(rho(t)=C*exp((a+I*b)*t 2 + (c+I*d)*t), ode [4]) ): 

# c is the shift of the pulse profile from extremum 

# of the modulation curve, d is the frequency shift of 

# the pulse from the center of the gain band 
expand(%/C/exp(t 2 *a)/exp(I*t 2 *b)/exp(t*c)/exp(I*t*d)): 

eql := collect (coeff(%,I), t ): 

eq2 := collect ( coeff(%%,I,0), t ): 
eq3 := coeff(eql, t 2 ): 

eq4 := coeff(eql, t): 
eq5 := coeff(eql, t, 0): 
eq6 := coeff(eq2, t 2 ); 
eq7 := coeff(eq2, t): 
eq8 := coeff(eq2, t, 0): 
allvalues( 
solve(subs(b=-a,{eq3=0,eq4=0,eq5=0,eq7=0,eq8=0}),{a,c,d,g,phi}) ); 



o2 o 

od ei ■= ( — p(t)) + gp(t) + S (— p(t)) +I<f>p(t) + IMt 2 p(t) 



eq6 :=Aa 2 -Ab 2 

{c = ~5, 4> = \ V2 VM, a = 1 V2 VM, g = ~ V2 Vm + 1 S 2 , d = 0}, 
{c=~6, .g = i\/2VM+i(5 2 , 4>=-^V2VM 7 a = -^v / 2v / M, d = 0} 
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We see, that there is not frequency shift of the pulse (d=0), but, as it was 
for amplitude modulation, the time delay appears ( c = — |) that changes 
the pulse duration as result of modulation detuning. The rise of the detuning 
prevents from the pulse generation due to saturated net-gain coefficient increase 



The pulse parameters behavior coincides with one for amplitude 
modulation 



5.3 Ultrashort pulse stability 

Now we shall investigate the ultrashort pulse stability against low perturba- 
tion C(t). The substitution of the perturbed steady-state solution in dynamical 
equation with subsequent linearization on £(£) results in 



— C(z, t) = a({z, t)+a p p(t) - l((z, t) + — 



— C(z, t) = a((z, t)+a p p(t) - K{z, t) + ^ C(z, t) - Mt 2 ((z, t), 



where a p is the perturbed saturated gain, which is obtained from the assumption 
about small contribution of perturbation to gain saturation process: 

>alpha[p] = op(2,convert( series( alpha[0]/(l + A + B), B=0,2), 
polynom)); 



a{) B 
a v = 



{1 + Af 



Here A= f_ a 2 dt, B= J aQdt and we neglected the high-order terms rela- 
tively perturbation amplitude. Note, that this is negative quantity. 

Let the dependence of perturbation on z has exponential form with increment 
A. Then 

>zeta(z,t) := upsilon(t)*exp(lambda*z):# perturbation 
ode[4] := expand( diff(zeta(z,t),z)/exp(lambda*z) ) = 
expand( (alpha*zeta(z,t) + alpha[p]*rho(t) - l*zeta(z,t) + 
diff(zeta(z,t),'$'(t,2)) - M*t 2 *zeta(z,t))/exp(lambda*z) ); 



ode 4 :=v(t)X = av(t)+ %^J- - lv(t) + (—^ v(t)) - Mt 2 v(t) 
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Now we introduce a v * where B*= f avdt, that allows to eliminate the 

V J— oo ' 

exponent from right-hand side of ode± (we shall eliminate the asterix below). 
This is equation for eigenvalues A and eigenfunctions v(t) of the perturbed laser 
operator. Stable generation of the ultrashort pulse corresponds to decaying of 
perturbations, i.e. A<0. For Gaussian pulse: 

>ode[4] := upsilon(t)*lambda = 

subs(rho(t)=sqrt(surd(M,4)*(alpha[0]-l-sqrt(M))/(sqrt(Pi)*(l+sqrt(M)))) 

*exp(-t 2 *sqrt(M)/2), 

alpha*upsilon(t)+alpha[p]*rho(t)-l*upsilon(t)+diff(upsilon(t),'$'(t,2)) 

-M*t 2 upsilon(t)); 

sol := dsolve( ode[4],upsilon(t) ); 



«fe 4 :=t,(t) A = a«(t)+ op j5HdW£Hz^£±£5) e (-V»Vfft-) 



-lv(t) + (^v(t))-Mt 2 v(t) 



J surd(M, 4)(l + VM-a ) ( _ 1/2 ^^) 

/ <A V" l + VM 

sot := v(t) ■■ 





(-A + a-Z- VM)tt( 1 /4) 

Ci WhittakerMf * ~ A + " " * , * , 

4 y/M 4' 


VMt 2 ) 


+- 


Ci2 WhittakcrW( 1 ~ X+ ^_~ l * 
4 VM 4' 


VMt 2 ) 



Vt 

As result (see above), we have 



v(t) = - 1 ^^ tt, + C HermiteH 2n (t vWM) e(-HF-) 

{<x-1-s/M-\)tt ( 4> 

for 

n= - ( 3 - 7^ && n is integer ' 



vTv7( Q0 -l-yTv7) e (- •"' ) 2 

u(t) = -^ ^L m + C HermiteH 2n+1 (t V7m) e ("MF) 

for 
n= " ( I _ W> && n is in teger, 
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>eql := n = - (l/4-(alpha-lambda-l)/(4*sqrt(M))): 
eq2 := n = - (3/4-(alpha-lambda-l)/(4*sqrt(M))): 
En := subs( 
C 2 =-surd(M,4)*(l+sqrt(M)-alpha[0])/(sqrt(Pi)*(l+sqrt(M))), 
C 2 *sqrt(Pi)/(M^/ 4 )) ):# ultrashort pulse energy (see above) 
lambda[l,n] = solve(eql, lambda);# increment 
lambda[2,n] = solve(eq2, lambda) ;# increment 



Ai n = -4 n \[~M - Vii? + a - I 



-AnVM -3VM + a-l 



The glance on the solution is evidence of absence of solution corresponding 
to Ai,o- For others modes we take into account, that a=l+ \[M. As result, 



all perturbation modes are unstable because of 



Ai,„ = -AnVM < (n =1, 2, ...), 



A 2 ,„ = -4nVM - 2VM < (n =0, 1, ...) 



It has to note, that for the amplified pulse 

Ai, „ = -4 n \[M - \[M + a - I, A 2 , „ = -4 n VM - 3 VM + a - I, 



Here we see the decrease of increment as result of n rise. There- 
fore, only "ground state" with A^o will be amplified predomi- 
nantly. This fact provides for Gaussian pulse generation 



(see, §). 

In the presence of detuning we have: 

>eql := n= - (l/4+(delta 2 -4*(alpha-lambda-l))/(16*sqrt(M))): 
eq2 := n= - (3/4+(delta 2 -4*(alpha-lambda-l))/(16*sqrt(M))): 

lambda[l,n] = simplify ( solve(eql, lambda) ); 
lambda[2,n] = simplify( solve(eq2, lambda) ); 
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Ai „ = -4 n VM - VM - - 5 2 + a - I 



\ 2n = -AnVM-'S^M- -S 2 + a-l 



Surprisingly, but, as it was above, our linear analysis predicts 
the pulse stability regardless of detuning S (because of a=l+ 



m + t) 



But for the pulse the detuning growth decreases the increment that does not 
favor the single pulse generation. 



5.4 Nonlinear processes: self-phase modulation and dy- 
namical gain saturation 

Among above considered effects only gain saturation by full pulse energy can 
be considered as nonlinear process, which, however, does not affect on the pulse 
envelope, but governs its energy. The time-dependent nonlinear effects, which 
can transform pulse profile, are self-phase modulation (SPM) and dynamical 
gain saturation. First one is the dependence of the field phase on its intensity 
and can play essential role in solid-state lasers ||. Second effect is caused by 
the change of the gain along pulse profile and is essential in lasers with large 
gain cross-sections and comparatively narrow gain band ||. 

At first, let analyze the presence of SPM, which can be considered as per- 
turbation of our master equation describing active phase modulation: 

>ode[5] := diff(rho(t),T(t,2)) + g*rho(t) + I*phi*rho(t) + 
I*M*t 2 *rho(t) _ i*beta*C 2 *exp(2*a*t 2 )*rho(t);# phase perturbation in the 
form -I*beta*C 2 *abs(rho(t)) 2 *rho(t), rho(t) has unperturbed profile 



ocfe 5 : = (^2 /»(*)) + 9p(t) +14>p(t) +IMt 2 p(t) -I(3C 2 e < 2 * 2 a > p(t) 



>expand( subs(rho(t)=C*exp((a+I*b)*t 2 ), 
ode[5])/exp(a*t 2 )/exp(I*t 2 *b)/C): 

series(%,t=0,3):# approximate solution as result of expansion 
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convert(%, polynom):# now we shall collect the coefficients oft and t 2 for 
real and imaginary parts 
eql := coeff(%,I): 

eq2 := coeff(%%,I,0): 
eq3 := coeff(eql,t 2 ); 
eq4 := coeff(eql,t,0); 

eq5 := coeff(eq2,t 2 ); 
eq6 := coeff(eq2,t,0); 
soil := solve(eq6=0, a);# solution for a 
sol2 := solve(subs(a=soll,eq3=0), b);# solution for b 
sol3 := solve(subs(b=sol2,eq4)=0, phi);# solution for phi 
# but NB that there is eq5 defining, in fact, g: 
solve( subs({a=soll, b=sol2},eq5)=0, g);# solution for g 



eq3 :=8ab + M -2aC 2 (3 



eq4 :=cl) + 2b-C 2 [3 



eq5 :=Aa 2 -Ab 2 



eq6 :— 2 a + g 



soil := q 

2 y 



sol2 := I M + 9C2 P 

4 9 



solS 



1 -M + gC 2 p 

2 9 



- C 2 13 + 1 VC 4 /3 2 + 8M, J C 2 13 - i VC 4 /3 2 +8M, 
-\° 2 P + l VC 4 f3 2 - 8 M, - j C 2 [3 - - ^C* f3 2 -8M 



Hence we have solution with perturbed parameters: 



P(t) 



II fJC 2 + ^p2 C 4 +8M , I (gC 2 +x/g 2 C 4 + 8M) ^ 2-, 

C e {{ 4 4 ' ' 
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i 1/3C 2 iJWd^+SM . 

and 5= — ^ 1 - j . As one can see, 



this solution has enlarged chirp and reduced pulse duration due 
to SPM 



Now we shall investigate the influence of dynamical gain saturation on the pulse 
characteristics in the case of active mode locking due to amplitude modulation. 
Let the contribution of the dynamical gain saturation can be considered as 
perturbation for Gaussian pulse (see above). The instant energy flux of such 
pulse is: 

>int(sol_0*exp(-t*(delta+sqrt(M)*t)/2),t): 
Energy := simplify( coeff( %, 

erf(l/2*sqrt(2)*M( 1 / 4 )*t+l/4*delta*sqrt(2)/(M( 1 / 4 ))))* 
(l+erf(l/2*sqrt(2)*M( 1 / 4 )*t+l/4*delta*sqrt(2)/(M( 1 / 4 ))) ) ); 



Energy := 



2 4 / + 4 VM + S 2 

The approximation of the small contribution of the gain saturation allows the 
expansion of energy into series on t up to second order: 

>series(Energy,t=0,2): 
convert (%, polynom): 
simplify (%); 

' (Al + iVM + S 2 -4a )e (_1/8 A ) 






Hence we have the modified gain coefficient a( 1 — lot), where a is the gain 
coefficient at pulse peak, 1$ is the pulse intensity for unperturbed solution. Note, 
that the additional term in brackets is resulted from the shift of pulse maximum. 
Then the master equation for perturbed solution (in our approximation!) is 

>ode[6] := diff(rho(t),T(t,2)) + g*rho(t) + delta*diff(rho(t),t) - 
(M*t 2 +epsilon*t)*rho(t);# here epsilon is alpha*I[0] 



G2 



o2 f\ 

ode 6 :=(-^p(t))+gp(t)+5(-p(t))-(Mt 2 + et)p(t) 



>dsolve(ode[6]=0, rho(t)); 



P(t) = 
M5 2 +t 
16 M&/*) ' 4' 4 M( 3 / 2 ) 



,1 4M 9 -M<5 2 +e 2 1 l(2Mt + e) 2 x , 1/9jm 
Ci WhittakerM(- 9 — — , -, - ^——>) e (- 1 /^) 



2Mt + s 



M( 3 / 4 ) 
MS 2 +t 
16 M^72) ' 4' 4 M( 3 / 2 ) 



,1 4Mg-MS 2 +e 2 1 l(2Mt + e) 2 , , 1/0 , t s 
. CI? WhittakerW(- 9 — Fn[m , -,- { „> )e^' 2 ^ 



2Mt + e 
M( 3 / 4 ) 

The comparison with above obtained result gives 



pUt) = C HermiteH 2n (t VVM + — ^ ) e { smvm 2 ) f or n= . ( 

2VM\/M 

-) && n is integer, 



16 M VM 



P2(t) = C HermiteH 2n +i(t V VM + — ^ ) e { sm^/m 2 ) f or n= . ( 

2\/MVM 

t - 4g r 6 + M 2 ^^ 2 ) &* - - ****■ 

In future we shall omit e 2 -terms. If we take the unperturbed intensity for 
calculation of perturbation action, the perturbed pulse energy is 

>C 2 *HermiteH(0, t*sqrt(sqrt(M))+epsilon/(2*sqrt(M*sqrt(M))) ) 2 * 
exp(-(4*M 2 *t 2 +4*M*t*epsilon)/(4*M*sqrt(M))-delta*t):# field 
# intensity profile for "ground state" 

(l/4-(4*g*M-M*delta 2 )/(16*M*sqrt(M))) = 
subs(epsilon=alpha*l[0] ,%) : 

int (% ,t=-infinity. .infinity) : 
En := simplify (%);# energy of the first "ground state" 



1/4 



(a Iq + S s/MY 



eV ' - (3/2) ) C 2 0F 

kn '■= 7TTT*. 

MW 4 ) 
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>l/4-(4*(alpha-l)*M - M*delta 2 )/(16*M*sqrt(M)): 

sol_a := solve(%=0, alpha): 
eq := sol_a - alpha[0]/(l+subs( alpha=%,En )) = 0:# energy balance con- 
dition 

sol_I := solve(eq, C 2 );# pulse peak intensity 



sol _I := 

M< 5 / 4 ) (4Z + 4 VM + S 2 - 4 oto) 



i laA ( 4/ m( 3 / 2 )+4J M l+In M S 2 +45 m( 3 / 2 ))' 

1 / 64 IPJ72) 



Pk (4 M( 3 / 2 ) + 4 Af I + M S 2 ) 



Now we plot the dependencies of the pulse intensity, width and maximum 
location versus detuning parameter S. 

>plot(subs( {alpha[0]=1.2,l=0.1,M=0.05},subs(I[0]=sol_0,sol_I) ), 
delta=-2..1.5, axes=BOXED, title='pulse intensity vs detuning'); 



pulse intensity vs detuning 




-0.5 

delta 



>eq := exp( -(4*M 2 *t 2 +4*M*t*epsilon)/(4*M*sqrt(M))-delta*t ) = 1/2: 
# we take into account pulse profile without epsilon 2 

sol := solve(eq, t): 

pulse_width := simplify( subs( epsilon 2 =0,sol[l] - sol[2] ) ): 

subs(epsilon=alpha*l[0] ,pulse_width) : 
subs({alpha=sol_a,l[0]=sol_0},%): # pulse width for perturbed solution 
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>plot(subs( {alpha[0]=1.2,l=0.1,M=0.05},% ), delta=-2..1.5, 
axes=BOXED, title='pulse width vs detuning'); 



pulse width vs detuning 




># pulse maximum location 
diff(exp(-(2*M*t+epsilon) 2 /(8*M*sqrt(M))-delta*t/2), t) = 0: 

solve(%, t); 



1 e+SVM 
~2 M 

>subs(epsilon=alpha*l[0],%): 

subs({alpha=sol_a,I[0]=sol_0},%): 
plot(subs( {alpha[0]=l. 2,1=0. 1,M=0.05},% ), delta=-2..1.5, 
axes=BOXED, title='pulse location vs detuning'); 
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pulse location vs detuning 




-2 - 



We see, that the main peculiarity here is the asymmetric depen- 
dence of the pulse parameters on 6. The pulse width minimum 
and intensity maximum don't coincide with 6=0 and the detun- 
ing characteristics have sharper behavior in negative domain of 
detuning 



Now, as it was made in previous subsection, we estimate the condition of the 
ultrashort pulse stability. Here we take into consideration e 2 -term, but this 
does not fail our analysis because of this term contribute only to pulse energy 
without shift pulse inside modulation window. For the sake of the simplification, 
we shall consider the contribution of the destabilizing field to dynamical gain 
saturation, but only in the form of the unperturbed peak intensity variation ( 
and perturbation of the saturated gain coefficient (see above) . Then the stability 
condition: 

>eql := epsilon = subs(I[0]=sol_0,sol_a*(I[0]+zeta)): 
# zeta is the peak intensity variation 

-(3/4-(4*(sol_a-lambda-l)*M+epsilon 2 -M*delta 2 )/(16*M*sqrt(M))): 
expand( solve(%=0, lambda) ) < 0;# perturbation increment from eigen- 
value problem 

subs(epsilon=rhs(eql),lhs(%)): 

animate3d(subs( {alpha[p]=0,alpha[0]=1.2,l=0.1},% ), 
delta=-2..1.5, M=0.005..0.05, zeta=-0.15..0.15, axes=BOXED, 
title='maximal increment of perturbation' 
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We can see, that the perturbation growth can destabilize the 
pulse as result of \6\ increase (compare with subsection "Ampli- 
tude modulation") 



Also, there is the possibility of ultrashort pulse destabilization near 6=0. So, 
the presence of dynamical gain saturation gives the behavior of the ultrashort 
pulse parameters and stability condition, which is close to the experimentally 
observed and numerically obtained (see, for example, |iq| . 
Now we try to investigate the influence of dynamical gain saturation in detail 
by using so-called aberrationless approximation Let the pulse profile keeps 



its form with accuracy up to n-order of time-series expansion, but n pulse 
parameters are modified as result of pulse propagation. Then the substitution 
of the expression for pulse profile in master equation with subsequent expansion 
in i-series produces the system of n ODE describing the evolution of pulse 
parameters. 



>fl := 

(z,t)— >rhoO(z)*exp(-a(z) 2 *t 2 +b(z)*t);# approximate pulse profile 
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f2 := (z,tau)- >rhoO(z)*exp(-a(z) 2 *tau 2 +b(z)*tau): 
master := diff(rho(z,t),z) = (alpha- l)*rho(z,t) - 
chi*alpha*rho(z,t)*int(rho(z,tau) 2 ,tau=-infinity..t) + 
diff(rho(z,t),t$2) - M*t 2 *rho(z,t) + delta*diff(rho(z,t),t);# approxi- 
mate master equation with dynamical gain saturation (parameter chi) , alpha is 
the gain coefficient before pulse front 

expand(lhs(subs({rho(z,t)=fl(z,t),rho(z,tau)=f2(z,tau) 
},master))*exp(a(z) 2 *t 2 )/exp(t*b(z)) - 
rhs(subs({rho(z,t)=fl(z,t),rho(z,tau)=f2(z,tau) 
},master))*exp(a(z) 2 *t 2 )/exp(t*b(z))):# substitution 
of the approximate solution 

convert( series(%,t=l/2*b(z)/(a(z) 2 ),3),polynom ):# expansion around 
peak at t=l/2*b(z)/(a(z) 2 ) 
eql := collect(%,t): 

eq2 := 
subs({diff(rhoO(z),z)=u,difT(a(z),z)=v,diff(b(z),z)=w },coeff(eql,t 2 )): 
cq3 := 

subs({diff(rhoO(z),z)=u,diff(a(z),z)=v,diff(b(z),z)=w },coeff(eql,t)): 
eq4 := 
subs({diff(rhoO(z),z)=u,difT(a(z),z)=v,diff(b(z),z)=w },coeff(eql,t,0)): 
sol := simplify (solve({eq2=0,eq3=0,eq4=0},{u,v,w}) ); 



fl := (z,t)^pO(z) e (-^) 2 * 2 +bW*) 



master :— — - p(z, t) 
oz 



(a - 


-l)p(z,t)- 


J — oo 



r\2 f) 

(qP P(*> *)) " Mt2 P(*. f ) + S (gl /»(*. f )) 
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sol :={u= -j P 0{z)( X a pO(zf V* e (1/2 ^ ) V2 

lim -crf(- — ))a(z) 

t->(-oo) 2 a(z) 

- 4 a a(z) 2 + 8 a(z) 4 + 4 / a(z) 2 - 2 p0(z) 2 b(z) x a e (1/2 ^ } 
-4b(z)(5a(z) 2 -4b(z) 2 a(z) 2 ) /a(z) 2 , 

w = -xap0(z) 2 e il/2 ^ l) - 2 c5a(z) 2 - 4b(z) a(z) 2 , 
1 -M + 4a(z) 4 
V ~~2 &(z) * 

Now try to find the steady-state points of ODE-system, which correspond to sta- 

tionary pulse parameters. For this aim let introduce substitution X = X e 2Mz)2 '■ 

>sol_main := 

{w = -chi*alpha*rho0(z) 2 -2*delta*a(z) 2 -4*b(z)*a(z) 2 , 

v = l/2*(M-4*a(z) 4 )/a(z), u = 

l/4*rhoO(z)*(-chi*alpha*rhoO(z) 2 *sqrt(Pi)*sqrt(2)*a(z)+4*alpha*a(z) 2 

-8*a(z) 4 -4*l*a(z) 2 +2*rho0(z) 2 *b(z)*chi*alpha+ 

4*b(z)*delta*a(z) 2 +4*b(z) 2 *a(z) 2 )/(a(z) 2 )}; 



sol _main := 

{u = - / oO(z)(-xapO(z) 2 y/TrV2&(z) +4aa(zf - 8a(z) 4 - 4Za(» 2 

+ 2 pO(z) 2 b(z) X a + 4 b(z) <5 a(z) 2 + 4 b(z) 2 a(z) 2 ) /a(z) 2 , 

w = -Xap0(z) 2 -2(5a(z) 2 -4b(z)a(z) 2 , v= 1 M ~ 4 ^) } 



>eql := subs( {rhoO(z) 2 =x,a(z)=sqrt(y),a(z) 2 =y,a(z) 4 =y 2 }, 
expand(numer( subs(sol_main,u) )/rhoO(z) )=0 );# here we eliminated the 
trivial solution 

eq2 := subs( a(z) 4 =y 2 ,numer( subs(sol_main,v) )=0 ); 

eq3 := subs( {rhoO(z) 2 =x,a(z) 2 =y},subs(sol_main,w)=0 ); 
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eql := 

-Xot-% \/tt v2 y/y + 4ay — 8 y — 4Zy + 

2 xb(z) xot + 4b(z) 5 y + 4b(zf y = 

eq2 := M -4y 2 = 
eq3 := -\ a x - 2 5 y - 4 b(z) j/ = 



>allvalues( solve({eq2,eq3},{y,b(z)}) 



1 x Q x + <5 VM 1 r— 

{h{z) = —2 7m '^=2^ } ' 

{b(z)= 2 V^ ^ = -2^ } 



Hence 



2 



6 

2 a 2 



8 
s/m 



, b = — | — ^ /It ( x i s the pulse intensity). So, the shift is 
that differs from the usual result (see above) as result 



X ■>' <v 
W 



of dynamical gain saturation (last term), which shifts the pulse maximum in 
negative side. This additional shift has obvious explanation. The gain at the 
pulse front is greater than one at pulse tail due to gain saturation. This shifts 
the pulse forward as hole. Pulse width is: 

>eq := exp( t*(-delta/2-chi*x*alpha/(2*sqrt(M))) - sqrt(M)*t 2 /2 ) = 1/2: 
sol := solve(eq, t): 

pulse_width := simplify ( sol[l] - sol[2] ); 



pulse width :■ 


^M 5 2 + 2 S VM X ax + X 2 a 2 x 2 + 8 ikf( 3 / 2 ) ln(2) 


M 



We can see, that there is the minimum of the pulse duration 
in negative domain of S that corresponds to result, which was 
obtained on the basis of perturbation theory 
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The pulse intensity pO 2 : 

> Intensity := simplify ( 

solve ( 

subs({y = l/2*sqrt(M),b(z) = 

-l/2*(delta*M+chi*sqrt(M)*x*alpha)/M},eql), x)[l]);# pulse intensity 



Intensity 



V^FM( 3 / 4 ) +Sy/M- J M (tt y/M + 2 V¥ MiWd + 4a - A^M - 41) 



X« 



So, we have the following dependencies for the pulse duration 

>subs(x=Intensity,pulse_width): 

animate(subs( {alpha=l. 2,1=0.1},% ), delta=-4..4, M=0.01..0.1, 
numpoints=200, axes=BOXED, title='pulse width vs detuning' 



25- 



20- 



15- 



pulse width vs detuning 




10- 



-1 1 

delta 



and pulse intensity 

> animate (subs ( {alpha=l. 2, 1=0. l,chi=l}, Intensity ), delta=-4..4, 
M=0.01..0.1, numpoints=200, axes=BOXED, view=0..0.45, 
title= 'pulse intensity vs detuning' 



71 



pulse intensity vs detuning 



0.4- 



0.3- 



0.2- 



0.1" 




The last step is the stability analysis. The stability of our solutions can be 
estimated from the eigenvalues of Jacobian of sol. 

>eql := subs({rhoO(z)=x,a(z)=y,b(z)=z},subs(sol_main,u)): 
eq2 := subs({rhoO(z)=x,a(z)=y,b(z)=z},subs(sol_main,v)): 
eq3 := subs({rhoO(z)=x,a(z)=y,b(z)=z},subs(sol_main,w)): 
m[l,l] := diff( eql,x): 

m[l,2] := diff( eql,y ): 
m[l,3] := diff( eql,z ): 
m[2,l] := diff( eq2,x ): 
m[2,2] := diff( eq2,y ): 
m[2,3] := diff( eq2,z ): 
m[3,l] := diff( eq3,x ): 
m[3,2] := diff( eq3,y ): 
m[3,3] := diff( eq3,z ): 
A := 
array([[m[l,l],m[l,2],m[l,3]],[m[2,l],m[2,2],m[2,3]],[m[3,l],m[3,2],m[3,3]]]):# 
Jacobian 

Now we find the eigenvalues A of Jacobian directly by calculation of determinant. 

>evalm(A-[[lambda,0,0],[0, lambda, 0], [0,0, lambda]]): 

numer( simplify(det(%)) ): 
sol := solve(%=0, lambda): 
>x := sqrt (Intensity): 

y := sqrt( sqrt(M)/2 ): 
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z := -delta/2-x 2 *alpha/(2*sqrt(M)): 
si := evalf( subs({chi=l,l=0.1,alpha=1.2},sol[l]) ): 
s2 := evalf( subs({chi=l,l=0.1,alpha=1.2},sol[2]) ): 
s3 := evalf( subs({chi=l,l=0.1,alpha=1.2},sol[3]) ): 

plot3d({Re(sl),Re(s2),Re(s3)}, delta=-4..2, M=0.01..0.1, 
axes=boxed,title='Re(lambda) for initial perturbation') 




One can see that the pulse with Gaussian-like form is stable in 
the region of its existence (see two previous figures) 



We have to note that the considered here perturbations belong to limited class 
therefore this criterion is necessary but not sufficient condition of pulse stability 
(compare with previous consideration on the basis of perturbation theory, where 
we analyzed not only pulse peak variation but also gain coefficient change). 



6 Nonlinear Schrodinger equation: construction 
of the soliton solution by means of the direct 
Hirota's method 



The Schrodinger equation is the well-known nonlinear equation describing the 
weak nonlinear waves, in the particular, the laser pulse propagation in fibers. 
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In the last case, a pulse can propagate without decaying over large distance due 
to balance between two factors: SPM and group delay dispersion (GDD). These 
pulses are named as optical solitons [jy]. The ultrashort pulse evolution obeys 
to the next master equation: 



d d 2 

/ (^P)= fc 2(^2P)+/? \p\ 2 P 



which is the consequence of eq_field from part 3 in the case of transition to 
local time t->t-z/c. The right-hand side terms describe GDD (with coefficient 
k 2 ) and SPM (with coefficient /3), respectively. 

It is very important to obtain the exact soliton solutions of nonlinear equa- 
tions. There are the inverse and direct methods to obtain such solutions. One 
of the direct methods is the so-called Hirota's method. The main steps of this 
method are: 1) the selection of the suitable substitution instead of the function 
p (see the master equation) , that allows to obtain the bilinear form of the evolu- 
tion equation; 2) the consideration of the formal series of perturbation theory for 
this bilinear equation. In the case of soliton solutions these series are terminated. 

The useful substitution for the nonlinear Schrodinger equation is 



p(z,t) = G(z,t)/F(z,t). 

Let suppose that F is the real function. It should be noted that we can make 
any assumption about p to satisfy the assumptions 1) and 2). Hirota proposed 
to introduce a new D-operator in following way: 



Dz m Dt n ab =((&)- (&))<"((& )-(&)) n a(M)b(*i,ti) 

After substitution of p in the terms of functions G and F we obtain two bilinear 
differential equations with regard to the new operator D: 

[I Dz + k 2 Dt 2 } G F = 
k 2 Dt 2 FF - [3 G G* = (1) 

The functions G and F can be expanded into the series of the formal parameter 

e: 
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G = e Gl + e 3 G3 + e 5 G5 ; F = 1 + e 2 F2 + e 4 F4 



Let substitute G and F into Eq. (1) and treat the terms with powers of e 
as independent, to get the infinite set of the differential equations relatively Gl, 
G3, ...; F2, F4, ... . These formal series are terminated only in the case when 
the master equation has exact 7V-soliton solution. For instance, the set of first 
six differential equations in our case is: 



I(f z Gl) + k 2 (§ p Gl)=0; 

2k 2 (-§^F2)-f3Gl G1* = 0; 

I (& GS) + k 2 (^ GS) + [IDz + k 2 Dt 2 ] G1F2=0: 

2k 2 (^F4) + Dt 2 F2F2 -f3{G3 Gl*+ Gl G3*) = ; 

1 (m G5 ) + fc 2 im G5 ) + [IDz + k 2 Dt 2 } (G3 F2 + G1F4) = 0; 

2 k 2 (^ F6) + k 2 Dt 2 (F4 F2 + F2 F4) - (3 (G5 Gl* + G3 G3* + Gl G5*) = 



For sake of the simplification of the very cumbersome manipulations we 
introduce the procedure for operator Dt m Dz n , which acts on the functions a 
and b. The lasts are the exponents (or linear combination of the exponents) in 
the form e v , where rj(z, t) is linear function. 

>restart: 

with (plots): 

6.1 Procedure Dt m Dz 11 

>Dt_Dz := proc (a,b,m,n) 

local Summa,k,r, result: 
Summa := 0: 

if (n>l) and (m<>0) then 
for k from 1 to n-1 do 

Summa := Summa+binomial(n,k)*(-l) n ~ fe+m * 

der( der(b,z,(n-k)),t,m)*der(a,z,k) + 

binomial(n,k)*(-l) n ~ fe *der(b,z,(n-k))* 

der( der(a,z,k),t,m) 
od: 

fi: 

if (n>l) and (m>l) then 

for r from 1 to (m-1) do 
for k from 1 to (n-1) do 
Summa := Summa+binomial(m,r)* 
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binomial(n,k)*(-l)"- fc + m - r * 
der( der(b,z,(n-k)),t,(m-r))* 
der( der(a,z,k),t,r); 
od: 
od: 
fi: 

if (m>l) and (nOO) then 

for r from 1 to (m-1) do 
Summa := Summa+binomial(m,r)*(-l) m ~ r+ ™* 
der( der(b,z,n),t,(m-r))*der(a,t,r) + 
binomial(m,r)*(-l) m-r *der(b,t,(m-r))* 
der( der(a,z,n),t,r); 
od; 
fi: 

if (m<>0) and (n<>0) then 
Summa := Summa+(-l) m+n *der(der(b,z,n),t,m)*a+ 

(-l) m *der(a,z,n)*der(b,t,m) + (-l)"*der(a,t,m)* 

der(b,z,n)+der(der(a,z,n),t,m)*b; 
fi: 

if (m=0) and (n>l) then 
Summa := Summa+(-l)™*der(b,z,n)*a+der(a,z,n)*b; 
for k from 1 to (n-1) do 
Summa := Summa+binomial(n,k)*(-l)"~ fe * 
der(b,z,(n-k))*der(a,z,k); 
od: 
fi: 

if (m=0) and (n=l) then 
Summa := der(a,z,l)*b-der(b,z,l)*a: 
fi: 

if (n=0) and (m>l) then 
Summa := Summa+(-l) m *der(b,t,m)*a+der(a,t,m)*b; 
for r from 1 to (m-1) do 
Summa := Summa+binomial(m,r)*(-l)" l ~ r * 
der(b,t,(m-r))*der(a,t,r); 
od: 
fi: 

if (n=0) and (m=l) then 
Summa := der(a,t,l)*b-der(b,t,l)*a: 
fi: 



76 



if (n=0) and (m=0) then 
Summa := a*b 
fi: 

result := combine(Summa,exp): 

end: 

The next procedure will be used for calculation of the derivative of e v (or 
combination of exponents) on t or z with further simplification of the obtained 
expression. 



6.2 Procedure der 



>der := proc (f,x,m) 

local difF,i,function: 

subs(etal=etal(x),etals=etals(x), 
eta2=eta2(x) ,eta2s=eta2s(x) ,f ) : 
difF := diff(%,x$m): 

if (x=t) then 

function := subs({diff(etal(x),x)=bl, 

diff(eta2(x),x)=b2,diff(etals(x),x)=bls, 
diff(eta2s(x),x)=b2s},difF) 

else 

function := subs({diff(etal(x),x)=al, 
diff(eta2(x),x)=a2,diff(etals(x),x)=als, 
diff (eta2s(x) ,x) =a2s} ,difF) 

fi; 

subs(etal(x)=etal,etals(x)=etals, 
eta2 (x) =eta2 ,eta2s (x) =eta2s , function) : 

if (m>l) then 

difF := subs({diff(al,x)=0,diff(a2,x)=0, 
diff(bl,x)=0,diff(b2,x)=0,diff(als,x)=0, 
diff(a2s,x)=0,diff(bls,x)=0,diff(b2s,x)=0},%) 
else 
combine(%) 
fi: 

collect(%,exp): 

end: 
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The next procedure is used to calculate an integral of e 11 (or combination of 
exponents) on t or z with further simplification of the expression. 



6.3 Procedure Integr 

>integr := proc (f,x,m) 

local intF,i,gl,gls,g2,g2s, function: 

intF := subs(etal=gl*x,etals=gls*x,eta2=g2*x, 
eta2s=g2s*x,f): 

for i from 1 to m do 

intF := int(intF,x); 
od: 

if (x=t) then 

x := t; gl := bl; gls := bis; g2 := b2; g2s := b2s; 
else 

x := z; gl := al; gls := als; g2 := a2; g2s := a2s; 
fi: 

intF; 

collect(%,exp): 

subs(bl*t=etal,bls*t=etals,b2*t=eta2, 
b2s*t=eta2s,al*z=etal,als*z=etals,a2*z=eta2, 

a2s*z=eta2s,%); 
end: 

Now, let try to obtain a first-order soliton for nonlinear Schrodinger equa- 
tion. 

>macro(Gs='G*',Gls='Gl*',G3s='G3*',G5s='G5*', as='a*',bs='b*',etaOs 
= 'etaO*'): 

Gl := exp(etal):# successful substitution! 
I*der(%,z,l)+k_2*der(%,t,2): #first from the equations set 
factor(%); 

e vl (Ial +k_2bl 2 ) 

>#as result we can find the parameter al 
al := I*k_2*bl 2 : 

als := -I*k 2*bls 2 : 
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>Gls := exp(etals):# conjugated to Gl 
GlGls := combine(Gl*Gls): 

F2 := beta/(2*k_2)*integr(%,t,2);#F2 from the second equation of 



the set 



■y Q e (r]l + etals) 

F2 := 



2 k_2{bl + bis) 2 

># But the next equation of the set results in 
I_Dz_Gl_F2 := I*factor(Dt_Dz(Gl,F2,0,l)): 

d_Dt2_Gl_F2 := k_2*factor(Dt_Dz(Gl,F2,2,0)): 
factor(I_Dz_Gl_F2+d_Dt2_Gl_F2); 
Dt_Dz(F2,F2,2,0); 



To obtain this result we use the trivial relationships: 

>etal=al*z+bl*t+etalO: 

eta2=a2*z+b2*t+eta20: 

Dz*exp(etal)*exp(eta2) = (al-a2)*exp(etal+eta2); 
Dt 2 *exp(etal)*exp(eta2) = (bl-b2) 2 *exp(etal+eta2); 



Dze^e^ 2 = {Ik_2bl 2 - a2) e (r > 1+r > 2 ) 



Dt 2 e" 1 e" 2 = (bl -b2) 2 e { 



As was shown above a= - i hi b 2 , hence the last term in third equation of set is 
equal to 0. So, we are to choose G3 = to satisfy third equation. Furthermore 
Dt 2 F2 from fourth equation of the set is ( ( 2 tb+bsY i fc g) 2 ) Dt 2 exp( 77+ r/_s). 
But in concordance with above obtained relationships this expression is equal 
to zero. So we can choose F4 = 0. Thus to satisfy other equations we can keep 
in the expansion of functions G and F only Gl, G3 and F2. So, the formal 
series are terminated. Since e is independent parameter we can take e=l. 
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>rho := Gl/(1+F2); 



p:-- 



e'l 1 



— e (v^+etals) 

k_2(bl +bls) 



>subs(etal=al*z+bl*t+etalO,etals=als*z+bls*t+etalOs, rho): 
soliton := expand(subs({bls=bl,beta=l,k_2=l/2},%));#the choice of k_2 
and beta is only normalization of the values in equation 

> soliton := 

exp(l/2*I*bl 2 *z)*exp(bl*t)*exp(etalO)/ 
(l+l/4*exp(bl*t) 2 *exp(etal0)*exp(etal0s)/(bl 2 )); 



p (l/2Ibl 2 z) p (bl t) pj,10 

soliton := - 



_ r e (bl t)\2 „r/10 e etalOs 

1 + 4 



bl 2 



># Now we check the obtained solution by substitution of one in dynamical 
equation 

I*der(rho,z,l)/rho+beta*exp(etal+etals)/((l+l/2*beta*exp(etal+etals)/ 
(k_2*(bl+bls) 2 ))) 2 +k_2*der(rho,t,2)/rho: 
simplify (%); 



All right! This is the exact solution of the Schrodinger equation. Physically bl 
has a sense of the inverse pulse duration. So it is real parameter. But what is 
the free parameter rjlO? Let r\10 is real and e vl ° = bl. Then 

>subs({exp(etal0)=2*bl, exp(etal0s)=2*bl},soliton): 
simplify (%); 



e (i/2 6i (izbi+2t)) bl 



l +e (2blt) 



The obtained solution is the so-called first-order soliton with duration -jj , 
amplitude hi and phase bl z/2: 



p(z, t) = bl ■ sech(bl ■ t)exp{i ■ bl 2 z/2) 



>plot(subs({z=0,bl=l},%),t=-5..5,axes=boxed, title= 'first-order soliton'); 



first-order soliton 




But what is about different values oirjlOl The choice of the different values 
of rjlO results in the obtaining of the so-called collapsing pulses, i. e. pulses 
with singularity in the dependence of their parameters on z. 

The described here procedure is available for the analysis of the higher-order 
solitons. For example, the substitution of Gl = e vl + e v2 causes the termination 
of the series on fifth equation (you can prove this statement by using of the 
described above procedures). The obtained solution depends on z and is called 
as the second order soliton. 



In the case of Schrodinger equation, the main feature of the described above 
method is the termination of the formal series for an arbitrary order of the so- 
lution. Such behavior results from a very rich mathematical structure of the 
dynamical equation: an existence of the infinitely many nontrivial symmetries 
and conservation laws, Painleve property and integrability by means of the in- 
verse scattering method. The non-decaying pulse-like solutions of the integrable 
nonlinear evolutional equations are called as the solitons . But as we shall see 



81 



later, some nonlinear equations have the soliton-like solutions, but do not be- 
long to integrable class. For example, for nonlinear Landau- Ginzburg equation 
(see next part) there is the soliton-like solution in the first order of the Hirota's 
method. But the second-order solution does not lead to the termination of the 
series. Such soliton-like solutions of the noninte grable dynamical equation are 
called as the 



quasi-solitons (or solitary waves) 



Nonlinear Landau- Ginzburg equation: quasi- 
soliton solution 



Here we shall consider a soliton-like pulse, which is generated in the continuous- 
wave solid-state laser due to power-dependent saturation of the diffraction loss 
in the presence of the field self- focusing in active medium |l^, [la]. The action 
of the saturable loss can be described by the real cubic nonlinear term. The 
energy dissipation due to spectral filtering can be introduced by means of the 
real second-order derivative on t. Then in the absence of SPM and GDD, the 
dynamical equation is a analog of the Schrodinger equation, but with pure real 
terms. 



d d 2 

— p(z, t) = gp(z, t) + t f 2 — p(z, t) + crp(z, tf 



Here g is the net-gain in the laser taking into account the gain and linear loss in 
the active medium, output loss, and diffraction loss. For sake of simplification, 
we shall suppose the normalization of time to the inverse bandwidth of the spec- 
tral filter tf (let £/= 2.5 fs, that corresponds to the full generation bandwidth of 
Ti: sapphire laser) and normalization of pulse intensity to the inverse intensity 
of the loss saturation a (the typical values of a are ~ 10~ 10 - 1CP 12 cm 2 /W). 
We will consider the steady-state pulse propagation, then -^ p(z, t) = 0. So, the 
master equation can be transformed to the well-known Duffing's type equation 
describing nonlinear oscillations without damping: 

>restart: 

with (plots): 

with(DEtools): 
ode := diff(rho(t),'$'(t,2)) + rho(t) 3 + g*rho(t); 



d 2 

ode := (— p(t)) + p(t) 3 + g p(t) 



82 



Its implicit solutions are: 

>sol := dsolve(ode=0,rho(t)); 



rp(t) i 

sol := / 2 d a-t- C2 = 0, 

J y / -2_a 4 -4g_a 2 +4_Cl ~ 

rp(t) i 

- 2 d a-t- C2 = 

^/-2 a 4 -4g a 2 + 4 C*i ~ 



The first integral of motion is: 

>numer(diff(lhs(sol[l]),t)): 

int_motion := simplify((op(l,%) 2 -op(2,%) 2 )/2); 



d 

int_motion := 2 (— p(t)f + p{t) A + 2g p(tf - 2 _ CI 



These equations describe the motion in the potential: 

>pot := simplify(op(2,int_motion)+op(3,int_motion)); 

pot:=p{tf + 2gp{t) 2 



The value of 4 *_ CI plays a role of the full energy of system. The dependence 
of potential on p for the different g is shown in the next figure: 

>plot3d(subs(rho(t)=rho,pot),g=0.05..-0.1,rho=-0.5..0.5,axes=boxed, 
title='Potential of pendulum'); 



S3 




We can see, that for g > there is one equilibrium state of pendulum in p = 
(stable state) , and for g < there are three one (unstable in p = and stable in 
p = +/- y/—g). Obviously, that in this system the different oscillation regimes 
is possible, that is illustrated by the phase portrait on the plane [y, z=d p/dt]. 



>sys := convertsys(ode = 0,[],rho(t),t,z); 

dfieldplot ( [diff (rho (t ) ,t) =z (t) ,diff (z (t ) ,t) =-rho(t) 3 -subs (g= 
[z(t),rho(t)],t=-2..2,rho=-0.5..0.5,z=-0.1..0.1, 
arrows=LARGE,axes=boxed,title='Nonlinear 
oscillations' ,color=black) ; 



-0.1,g)* rho(t)], 



sys := [[YP l = z 2 , YP 2 = - Zl 3 - g Zi ] 



[zi = p(t), Z2 = tt p(t)], undefined, 
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Nonlinear oscillations 
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The quasi-soliton solution of the initial equation corresponds to the oscilla- 
tion around the stable equilibrium state with infinite period. In this case the 
full energy is equal to 0. Then _ CI = and the motion begins from p = at 
t = - oo: 

>plot({sqrt(subs({rho(t)=rho,g=-0.1},-pot))/2,-sqrt(subs( 
{rho(t)=rho,g=-0.1},-pot))/2}, 

rho=0..0.45,axes=boxed,labels=['rho(t)','drho(t)/dt'],color=red, 
numpoints=200); 
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rho(t) 
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The amplitude of quasi-soliton can be found from the integral of motion (i. 
e. pulse maximum correspond to d p/dt = and _C1 = in the integral of 
motion). 

>rho0 := solve(factor(pot)/rho(t) 2 =0,rho(t)); 



pi) := ypTg, -^2< 



The explicit integration of the solution produces: 

>assume(g<0): 

soll_a := mimer(value(subs(_Cl=0,lhs(sol[l])))); 
soll_b := numer(value(subs(_Cl=0,lhs(sol[2])))); 



soil _a := p(t) y/-2 p(t) 2 - 4 g~ arctanh(2 



V=FV-2pW 2 - 4 5 



=)- 



t y/-2p(t) 4 - 4 5 >(t) 2 V^F - _C2 ^-2p(ty - 4 ff >(t) 2 v^ 



soll_b := -p(t) v /-2p(t) 2 -4 5 ~arctanh(2 ^=== = ] 

- t V-2 p(t) 4 - 4 <f p(t) 2 ^f - _ C2 V-2 p(t) 4 - 4 «f p(i) 2 v^7 
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Make some transformations: 

>sol2_a := 

numer(simplify(expand(soll_a/op(l, soil _a)), radical, symbolic)); 

sol2_b := numer(simplify(expand(soll_b/op(l, soil _b)), radical, symbolic)); 



sol2_a := arctanh(2 — ) - ty/-g~ - _C2 -J^i 



sol2_b := arctanh(2 =) + ty/-g~ + _C2 \f- 

y^gy/-2p(t) 2 -4:g~ 



When at the pulse maximum p(0) = p m ax, the value of _ C2 can be found as: 

>i_C2_a := solve(subs({t=0,rho(t)=rho_max},sol2_a)=0,_C2); 
i_C2_b := solve(subs({t=0,rho(t)=rho_max},sol2_b)=0,_C2); 



arctanh(2 ) 

■J— g~ \ —2 rho max — 4g~ 

i C2 a:= ■ — V / ~ 

V-9 

n~ 

arctanh(2 ■ 



J — g~ \ — 2 rho max -4o" 

i_C2_b := ■ — V / ~ 

V-9 



Let tanh(i_C2_a*sqrt(-g)) = tanh(-i_C2_b*sqrt(-g) ) = -/+ v, then from an 
expression for the tangents of sum of arguments: tanh(a + b) = tanh(a) + 
tanh(b) / '(1 + tanh(a)tanh(b)), the equations can be transformed as: 

>sol3_a := 

simplify (tanh(op(l,sol2_a))) = (tanh(-op(2,sol2_a))+upsilon)/ 

(l+tanh(-op(2,sol2_a))*upsilon); 

sol3_b := 
simplify (tanh(op(l,sol2_b))) = (tanh(-op(2,sol2_b))+upsilon)/ 
(l+tanh(-op(2,sol2_b))*upsilon); 



,„ o 9~ t&nh{tV^T) + 

sou a := 2 — 



V^iT v / -2/oW 2 -4.9~ 1 + tanh(i V^T) v 
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sol3 b:=2 



g~ — tanh(£ \/ —g~) + v 

VFT \Z-"2p(t) 2 -4.g~ ~ 1 - tanh(t V^F) f 



>sol4_a := solve(sol3_a, rho(t)); 
sol4_b := solve(sol3_b, rho(t)); 



S0I4 _a := 2 M ^ ; — ; ^- — ; , —2 



%l-l + i>%l + i; 



%1 - 1 + v %1 + v 

%l:=(e( 4 ^?-))2 



sol4_b := 2 



V2 5 ~ -2<fi; 2 e(*v^r) J2g~-2g t-- , 



-,,2 pN-s") 



-%1 + 1 + u %1 + u 



-%1 + 1 + v %1 + u 



Now, we must note, that the transit p max -> pO (see above) corresponds to 
the transit v > oo. Then the final solutions result from the next operations: 

>sol_fin_l := limit(sol4_a[l],upsilon=infinity); 
sol_fin_2 := limit(sol4_a[2],upsilon=infinity); 

sol_fin_l := limit(sol4_b[l],upsilon=infinity); 
sol_fin_2 := limit(sol4_b[2],upsilon=infinity); 



sol fin 1 : 





_ n y/-2g-e^-«> 



sol_ftn_2 := -2 






sol_fin_l := 2 



e (2tV=9l +~T 



sol_fin_2 := -2 






The positive roots satisfy to the initial condition and the result is the quasi- 
soliton pulse with sech - shape envelope, which has the duration 



-<J 



and 



the peak amplitude pO = \f—2g . The pulse intensity profile is shown in the 
next figure ( a = 10 -11 cm 2 /W, I is the time normalized to tf): 



>animate(evalf(subs(t=iota/(2.5e-15),sol_fin_l) 2 *100), 
iota=-le-13..1e-13,g=-0.05..-0.01,frames=50, 
axes=boxed,color=red,labels=['time, fs','rho 2 , GW/cm 2 '], 
title='Pulse envelope'); 



Pulse envelope 




-8e-14 



-4e-14 



2e-14 
time, fs 



6e-14 



1e-13 



So, there is the close analogy between an ultrashort pulse gen- 
eration in the laser with power-dependent loss saturation and 
an oscillation of nonlinear pendulum. The considered solution 
has the character of steady-state quasi-soliton 



In the next part we are going to demonstrate some analogies of the high-order 
soliton behavior (laser breezers). 



8 Autooscillations of quasi-solitons in the laser 



The main efforts in the design of the modern femtosecond lasers aim to the 
stabilization of the ultrashort pulse parameters. For example, as it was shown 
[[L4J, the pulse destabilization are very important factors limiting the operation 
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of laser with so-called fast saturable absorber (the bleaching of diffraction loss 
due to Kerr-lensing is the example of such saturable absorber, see previous 
part). The destabilization of the laser quasi-soliton can produce its destruction 
or regular/nonregular autooscillations. The last can be considered as the analog 
of a high-order soliton formation (see part VI), which is reminiscent of the 
breezers of nonlinear dynamical equation. We shall consider the master laser 
equation, which joins the Landau - Ginzburg and Schrodinger equations. The 
main nonlinear factors now are SPM and power-dependent loss saturation. GDD 
and spectral filtering will be taken into consideration, too. 

>restart: 

with(DEtools): 
with (plots): 
master_l := diff(rho(z,t),z) = 

alpha*rho(z,t)-gamma*rho(z,t)+I*phi*rho(z,t)+tf 2 *diff(rho(z,t),'$'(t, 
2))+I*k_2*diff(rho(z,t),'$'(t,2))+sigma*rho(z,t) 2 *conjugate(rho(z,t)) 
-I*beta*rho(z,t) 2 *conjugate(rho(z,t)); 



master 1 := — — p(z, t) = 
~ az 



£)2 pil 

ap{z, t) -ip(z, t) + I(f>p(z, t) + tf (— p{z, t)) + Ik_2 (^ p(z, *)) 



+ ap{z, ty {p(z, t)) - I0p(z, ty {p(z, t)) 



Here </) is the phase delay on the full cavity round trip, a and 7 are the gain 
and loss coefficients, respectively. The general exact solution of this equation is 
not known, but there is the quasi-soliton solution in the following form: 

>fl := (z,t)- >rhoO(z)*sech(t*tau(z)) 1+/ *P s ^); 



fl := {z, t) -» pO(z)sech(tT(z)) il+I ^ z)) 



Here pO is the pulse amplitude, r is the inverse pulse width, V is the chirp. 
This solution obeys the condition of steady-state propagation, when J^ p = 0, 
i.e. the pulse parameters are constant. To describe the nonstationary pulse 
propagation, we shall use the aberrationless approximation: the changes of the 
pulse parameters do not cause the aberration of the pulse form. Next step is 
the substitution of fl into master _1 with following expansion in t - series. The 
coefficients of the expansion produce the set of ODE for the evaluating pulse 
parameters. 
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>assume(tau(z) ,real) : 
assume(t,real): 

expand(lhs(subs(rho(z,t)=fl(z,t),master_l))-rhs(subs(rho(z,t)=fl(z,t), 
master_l))): 

eq := 
subs( 
{alpha=alpha(z) ,diff (rhoO(z) ,z) =a,diff(tau(z) ,z) =b,diff(psi(z) ,z) =c} , 
convert (series(%,t=0, 3) ,polynom)): 
assume(rhoO(z) ,real) : 

eql := evalc(coeff(eq,t 2 )): 
eq2 := evalc(coeff(eq,t)): 
eq3 := evalc(coeff(eq,t,0)): 
eq4 := coeff(eql,I): 

eq5 := coeff(eql,l,0): 

eq6 := coeff(eq3,I): 
eq7 := coeff(eq3,I,0): 
solve({eq4=0,eq5=0,eq6=0,eq7=0},{a,b,c,phi}): 

sys := 
diff(rhoO(z),z)=subs(%,a),diff(tau(z),z)=subs( 
%,b),diff(psi(z),z)=subs(%,c): 

The obtained system sys have to be supplemented with equation for gain evo- 
lution: 

>sys := 

{ 

%,diff(alpha(z),z)=alpha(z)*exp(-2*xi*rhoO(z) 2 /tau(z)-l/Tr-Pump) + 
Pump*alphamx*(l-exp(-l/Tr-Pump))/(Pump+l/Tr)-alpha(z)}; 
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d 
sys ■= {-Q-: pO(z ) = 

a(z~)pO(z~) + <j p0(z~) 3 -jpO(z~) - 
tf pO{z~) t(z~) 2 + k_2pO(z~) i/>(z~) t{z~) 2 , 

d ( ~\ 
oz 

a P 0{z~f t(z~) - 2 tf 2 t{z~) 3 + 3 k_2 r{z~f ijj(z~) + tf 2 ijj(z~) 2 t(z~) 3 , 

oz 
-2 tf tp(z~) t(z~) 2 -4k_2tp(z~) 2 T(z~) 2 -2(3p0(z~) 2 -Ak_2r{z~) 2 

-2iP(z~)ap0(z~) 2 -2tf 2 ib(z~) 3 T(z~) 2 , 

, _. ( _ 2 epo(Q 2 i _ P ) Pumpalphamx(l-e(-Ty~ Purn Py) 
a(z )e v T < 2 ) Tt fi H = a(z )} 

Pwmp + — 
It 

Here Tr is the gain relaxation time normalized to the cavity period, Pump is 
the dimensionless pump (see part 9), £ is the inverse gain saturation energy. 

At first let's find the parameters of a steady-state quasi-soliton, viz. the 
solution of sys independent on z. 

>f[l] := 
-2*tf 2 *psi*Tau - 4*k_2*psi 2 *Tau - 4*k_2*Tau - 2*tf 2 *psi 3 *Tau - 2*sigma*Phi 
*psi - 2*beta*Phi;# Tau = tau 2 , Phi = rho 2 

>f[2] := alpha + sigma*Phi - gamma - tf 2 *Tau + k_2*psi*Tau; 

>f[3] := 3*k_2*Tau*psi + tf 2 *psi 2 *Tau - 2*tf 2 *Tau + sigma*Phi; 

/i := -2 tf 2 i>T - 4 k_2 f T - 4 k_2 T - 2 tf tp 3 T - 2 a $ ip - 2 /3$ 



h ■= a + a $ - 7 - tf T + k_2 i\) T 



h :=3k_2ipT+tf 2 fT-2tf 2 T + a<P 

>soll := solve({f[2]=0,f[3]=0},{Tau,Phi}); 
>subs({Tau=subs(soll,Tau),Phi=subs(soll,Phi)},f[l]): 
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> simplify (%): 

>numer(%)/2/(-alpha+gamma): 
>so!2 := solve(%=0,psi); 



soil := {$ = _("-7)(3t_^ + t/ 2 ^-2 f / a ) i T 



a — 7 



(2k_2ip + tf Z Tjj 2 -tf z )a 2k_2ip + tf Z TP 2 -tf 



2 ,1,2 



72} 



sol2 



1 


-3/3fc_ 


_^ + 3t/ 2 cr 


+ y / 9/3 2 fc_ 


_2 2 


-2(3k_ 


_2tf 2 o 


+ 9 tf 4 a 2 


+ 8/3 2 


tf 4 


+ 8k_ 


_2 2 


a 2 


2 


-3/3fc_ 


_2 + 3tf 2 a 






(3tf 2 + k_2a 














1 


-^9(3 2 k_ 


2 2 


-2(3k_ 


_2tf 2 o 


+ 9 tf 4 a 2 


+ 8/3 2 


tf 4 


+ 8k_ 


2 2 


a 2 



2 (3tf 2 + k_2a 

Normalization of the time to tf and the intensity to j3 allows simple plotting 
of the obtained result: 



>plot({subs( {beta=l,tf=l,sigma=l},sol2[l] ), 
subs( {beta=l,tf=l,sigma=l},sol2[2])}, 
k_2=-100..10,axes=boxed,view=-10..1, title= 'chirp vs. GDD' 

chirp vs. GDP 




-80 



-60 -40 
k 2 



20 



To plot the pulse duration and its intensity we have to take into account the 
dependence of a on the pulse energy. However for the sake of simplification we 
suppose a=const< 7 (see previous section). Then 
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>plot (subs( {beta=l ,tf= 1 ,sigma=l ,gamma=0.05, 
alpha=0.04},subs(psi=sol2[2],l/sqrt(subs(soll,Tau)))), 
k_2=-100..10,axes=boxed, title='pulse width vs. GDD'); 



pulse width vs. GDP 



100 




-1 00 -80 



-60 -40 
k 2 



Only solution with ip>Q in the region of the negative GDD and -0<O for the 
positive GDD (green curve) has a physical meaning because it corresponds to 
the positive square of the pulse width 1/ r. 



The existence of the pulse duration minimum in the vicinity of 
zero GDD pays a very important role in the pulse shortening 
technique based on the creation of appropriate negative GDD 
in the laser cavity. 



Now, we shall consider the evolution of ultrashort pulse parameters on the 
basis of the obtained system of ODE. We suppose to solve the system by the 
standard operator DEplot. Let normalize a and (3 to 1,7* 10~ 12 cm 2 /W, times 
to tf (2,5 fs for Ti: sapphire laser), then £ = 0.0018. The fundamental step is 



the assumption about saturation of the Kerr nonlinearity: a = o"o(l 



Pa a o 



), 



(3 = /3o(l - Po 2 ), where Co and j3q are the unsaturated nonlinear parameters. 

>#procedure for solution of the obtained system of ODE 

ODE_plot := proc(alphamx, gam, sigma0,beta0,Tr, Pump, xi,disp,tg,n) 
sigma := sigmaO*(l - sigma0*rho0(z) 2 /2 ): 
beta := beta0*(l - beta0*rho0(z) 2 /2 ): 
sys := [D (alpha) (z) = 
alpha(z)*exp(-2*xi*rhoO(z) 2 /tau(z)-l/Tr-Pump) + 
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Pump*alphamx*(l-exp(-l/Tr-Pump))/(Pump+l/Tr)-alpha(z), 

D(psi)(z) = 
-4*disp*tau(z) 2 -2*tg 2 *psi(z)*tau(z) 2 -4*disp*psi(z) 2 *tau(z) 2 - 
2*beta*rho0(z) 2 -2*beta*rho0(z) 2 *psi(z) 2 -2*tg 2 *psi(z) 3 *tau(z) 2 , 

D(rhoO)(z) = 
sigma*rho0(z) 3 -gam*rho0(z)-tg 2 *rho0(z)*tau(z) 2 +alpha(z)*rho0(z)- 
disp*rhoO(z)*psi(z)*tau(z) 2 , 

D(tau)(z) = 
-2*tg 2 *tau(z) 3 +sigma*rho0(z) 2 *tau(z)+3*disp*tau(z) 3 *psi(z) + 
beta*rho0(z) 2 *psi(z)*tau(z)+tg 2 *psi(z) 2 *tau(z) 3 ]: 

DEplot(sys,[rhoO(z),tau(z),psi(z),alpha(z)],z=0..n,[[rhoO(0)=0.001, 

tau(0)=0.01,alpha(0)=0,psi(0)=0]],stepsize=l,scene=[z,rhoO(z)], 

method=classical[foreuler],axes=FRAME,linecolor=BLACK): 



end: 



The next figure demonstrates the autooscillations of pulse amplitude (quasi- 
breezer behavior). 

>display(ODE_plot(0.5,0.05, 10,1,300,0.0004,0.0018,-10,1, 15000)); 
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The character of the pulse evolution strongly depends on the correlations be- 
tween system's parameters. For example, next figure demonstrates the pulse 
stabilization due to negative dispersion growth. 



95 



>display(ODE_plot(0.5,0.05, 10,1,300,0.0004,0.0018,-20,1, 15000)); 



0.025 

0.02 

fjO.015 

0.01 

0.005 



o 
o 




2000 4000 6000 8000100001200014000 



But the pumping growth produces the irregular autooscillations. 
>ODE_plot(0. 5,0.05, 10, 1,300,0.00047,0.0018,-20,1,5000)); 
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The next parameter's set produces the chaotic oscillations: 

>display(ODE_plot(0.5,0.05, 1,0,300,0.0004,0.0018,0, 1,8000)); 
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The pulse parameters (amplitude and inverse pulse duration) corresponding 
to irregular oscillations can be shown on the basis of the iteration procedure, 
which realizes the direct Euler method for the solution of the described above 
ODE system. 

># Attention! This block can take a lot of CPU-time 
iterations := proc(iter) 
alphamx := 0.5: 
gam := 0.05: 
sigmaO := 1: 

betaO := 0: 
Tr := 300: 
Pump := 0.0004: 
xi := 0.0018: 
disp := 0: 
tg := 1: 

rho0n:=0.001; 
taun:=0.01; 
alphan:=0; 
psin:=0; 

for m from 1 to iter do 
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rhoOold:=rhoOn; 

tauold:=taun; 
psiold:=psin; 
alphaold : = alphan ; 
sigma := evalhf(sigmaO*(l - sigma0*rho0old 2 /2 )): 
beta := evalhf(betaO*(l - beta0*rho0old 2 /2 )): 
alphan := evalhf( 
alphaold*exp(-2*xi*rho0old 2 /tauold-l/Tr-Pump) + 
Pump*alphamx*(l-exp(-l/Tr-Pump))/(Pump+l/Tr)); 

psin := 
evalhf(psiold-4*disp*tauold 2 -2*tg 2 *psiold*tauold 2 - 
4*disp*psiold 2 *tauold 2 -2*beta*rho0old 2 - 
2*beta*rho0old 2 *psiold 2 -2*tg 2 *psiold 3 *tauold 2 ); 

rhoOn := 
evalhf(rho0old+sigma*rho0old 3 -gam*rho0old-tg 2 *rho0old*tauold 2 
+alphaold*rho0old+disp*rho0old*psiold*tauold 2 ); 

taun := 
evalhf(tauold-2*tg 2 *tauold 3 +sigma*rho0old 2 *tauold+3*disp*tauold 3 * 
psiold+beta*rho0old 2 *psiold*tauold+tg 2 *psiold 2 *tauold 3 ); 

if m = iter then pts := [rhoOn, taun] fi; 

od; 

pts 
end: 

PLOT(seq(POINTS(iterations(i)), i=9800 .. 10000), SYMBOL(POINT)); 
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This is the so-called strange attractor, i. e. the chaotic attracting manifold: the 
pulse parameters are changed chaotically but within the limited range. 

Thus, the pulse osc illations in the Kerr -lens mode-locked laser were analyzed 
(see Jl5J, and arXiv: \physics/000902(\i . 



The oscillations accompany the negative dispersion decrease and 
the pump growth and close connect with the nonlinear factors in 
the system. The regular oscillation is the analogue of the high- 
order soliton propagation and the irregular one is the analogue 
of the nonlinear breezer 



9 Numerical approaches: ultrashort pulse spec- 
trum, stability and multipulsing 

9.1 Kerr-lens mode- locked Cr:LiSGaF-laser with the Ra- 
man self-scattering in active medium 

The aberrationless approximation demonstrates the stability loss in the vicinity 
of zero dispersions. It is necessary to interpret the phenomenon. Moreover, we 
did not consider a lot of the lasing factors, which affect the ultrashort pulse 
dynamics in sub-100 fs region, viz. higher-order dispersion, stimulated Raman 
self-scattering (see Part II), strong fast absorber saturation etc. It is clear, 
that in order to take into account these phenomena we need the numerical 
simulations beyond the computer algebra abilities. However, Maple can help 
in the preparation of the source code and in the interpretation of the obtained 
results. 

Now we describe the simplest generation model, which is highly useful for the 
numerical simulations. This model is based on the generalized Landau- Ginzburg 
equation: 

>restart: 

>with(codegen,fortran): 
>with('linalg'): 
>with(stats): 
>with(plots): 

>masterl := diff(a(z,t),z) = 
(alpha(z) - rho - gamma/(l+sigma*Phi(z,t)))*a(z,t) + 
diff(a(z,t),t$2) + sum('(-I) fe+1 *D[k]*diff(a(z,t),t$k)','k'=2..N) - I*Phi(z,t)*a(z,t); 
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# Phi is the field intensity, gamma is the fast absorber modulation depth, 
sigma is the inverse saturation intensity, we used the normalization of the field to 
the self-phase modulation coefficient and the time to the inverse gain bandwidth, 
D[k] are the dispersion coefficients absorbed the (f /factorial(k)) factors 



masterl :— -^ &(z, t) = (ct(z) — p — - — ; — f — — ) a(z, t) + ( JW a(z, £)) 



l + a<f>(z, t)- 



N 



^(_/)(fc+i)D fe diff(aO, t), t$k)\ -I$(z, t)&(z, t) 



1.—2 



For the gain evolution we have: 

>master2 := diff(alpha(z),z) = P*(a - alpha(z)) -alpha(z)*Int(Phi(z,t),t= 
T[cav]/2..T[cav]/2)/E[s] - T[cav]*alpha(z)/T[r]; 



master2 :- 


oz 


= P(a- 


-a(z))- 


a(z) 


rl/2T cav 

1 $(z,t)dt 

-1/2T„ 


•*■ cav OtyZJ 

T r 




E s 



Here a is the maximal gain for the full population inversion, P is the di- 
mensionless pump (Pump from the previous section), T cav and T r are the cavity 
period and the gain relaxation time, respectively, E 8 is the gain saturation en- 
ergy taking into account the accepted normalization of time and intensity. The 
best methods for the solution of the system (masterl, master 2) is the split-step 
Fourier method. Then in the Fourier domain we escape to calculate the partial 
derivatives: 

>op(2,rhs(masterl)); 
>inttrans[fourier](%, t, omega) *F (omega); 

>print('F(w) takes into account the shape factors for the gainband and 
output coupler spectral profiles'); 

>subs( N=6,op(3,rhs(masterl)) );# we confine the maximal disper- 
sion order 

>inttrans[fourier](evalc( % ), t, omega): 
> factor (%); 

d 2 



-co 2 fourier(a(z, t), t, uj)F(uj) 
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^(-/)( fe+1 >D fe diff(a(z, t),t$k) 



fc=2 



-Iuj 2 fourier(a(z, t), t, w) (D 6 w 4 + D 5 lu 3 + D 4 uj 2 + D 3 cj + D 2 ) 



To define the parameters of the simulation we use the experimental data for 
Cr:LiSGaF-laser kindly presented by Dr. I.T. Sorokina and Dr. E. Sorokin and 
written in the corresponding *.txt-files. NOTE! ALL DATA FILES HAVE TO 
BE PLACED IN YOUR CURRENT DIRECTORY: 



>currentdir(); 



"G:\\Maple 6" 



We used the next experimental setup, which is typical for the Kerr-lens 
mode-locked lasers (high reflective plane and spherical (focal length f=5 cm) 
mirrors HR, chirped mirrors Ch, output mirror OC): 



HRorCM 3 CriiSGaF 



CM, 




The first step is the analysis of the gain band, which allows to define the gain 
bandwidth and the parameters of the time normalization (tf from the previous 
section) and the frequency normalization (this is the gain bandwidth). 



101 



>gain_x := readdata('gain_x.dat',l,float):# wavelength 
>gain_y := readdata('gain_y.dat',l,float):# gain cross-section 
>n := vectdim(gain_y): 

>p\ot([[gain_x[k], gain _y[k]] $k=l..n], 
color=red,view=0..3e-20,axes=boxed,title='gain band profile'); 

>max_cross_section := max (gain _y[k] $k=l..n);#gain band max- 
imum 

>half_cross_section := evalf(max_cross_section/2);#half of the gain 
band maximum 

>P := array(L.n): 
>Q := array(l..n): 

>for i from 2 to n do 

>P[i] := evalf(2*Pi*3*lelO/(gain_x[i]*le-7)):#transition from wavelength 
to frequency 

>Q[i] := [P[i],gain_y[i]]: 

>if gain_y[i]=max_cross_section then X_max := P[i] else fi: 
>if gain_y[i]>half_cross_section and gain_y[i-l]<half_cross_section 
then X_half_l := evalf((P[i-l]+P[i])/2) else fi: 

>if gain _y[i]<half_ cross _ section and gain_y[i-l]>half_cross_section then 
X_half_2 := evalf((P[i-l]+P[i])/2) else fi: 
>od: >print('position of the gain maximum:'); 

>X_max; 

>print('position of the first half of maximum:'); 

>X_half_l; 

>print('position of the second half of maximum:'); 

>X_half_2; 

>print ('bandwidth:'); 

>bandwidth := evalf(abs(X_half_2-X_half_l)); 

> print ('inverse bandwidth for Gaussian approximation [s]:'); 

>minimal_pulse_width:=evalf(4*ln(2) /bandwidth); 

>for i from 1 to n do 

>P[i] := evalf((2*Pi*3*lelO/(gain_x[i]*le-7)-X_max)/bandwidth): 

>Q[i] := [P[i],gain_y[i]]: 

>od: 

>j := T: 

>plot([Q[j] $j=l..n],axes=BOXED, color=red, title='gain cross section ver- 
sus normalized frequency', view=0..3e-20,color=blue); 
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gain band profile 
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max cross section :— .2926 10 
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half _cross _section := .1463000000 10" 19 



position of the gain maximum 



.224590912210 



id 



position of the first half of maximum 



.2476565668 10 16 



position of the second half of maximum 



.1994706028 10 16 
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bandwidth 



bandwidth := .481859640 10 15 



inverse bandwidth for Gaussian approximation [s] 
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minimal _pulse _width := .575393432510 
gain cross section versus normalized frequency 
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The similar manipulations are performed for the output coupler. The fit- 
approximation gives the spectral shape factor for the output mirror. 

>Out_x := readdata('out_x.dat',l, float): 
>Out_y := readdata('out_y.dat',l,float): 
>n := vectdim(Out_y): 

>plot([[Oui_x[/c],Oui_j/[fc]]$k=l..n],color=red,axes=boxed, 
title='experimental reflection profile'); 
>P0 := array(l..n): 
>Q0 := array(l..n): 
>for i from 1 to n do 
>P0[i] := 
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evalf( (2*Pi*3*lelO/ (Out _x[i]*le-7)-X_max) /bandwidth) :#transition from wave- 
length to frequency 

>Q0[i] := evalf(l-Out_y[i]): 

>od: 

>fl := plot([[P0[fc], Q0[fc]] $k=l..n],color=red):#experimental output loss 
profile in the dimensionless frequency domain 

>eqO := fit[leastsquare[[x,y], 
y=a*x 6 + b*x 5 + c*x 4 + d*x 3 + e*x 2 + f*x + g]]([[P0[fc]$A; = l..n], [Q0[k]$k = 
l..ra]]);# fit-approximation 

>Q0 := [subs(x=PO[k],rhs(eqO)) $k=l..n]: 

>f2 := plot([[P0[fc],Q0[fc]] $k=l..n],color=blue): 

>display(fl,f2,axes=boxed,title='comparison of the experimental and ap- 
proximated profiles'); 

experimental reflection profile 
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eqO :=y = 1.800575680 x 6 - .05353754671a; 5 - .1603422410 a; 4 - .002641254098 a; 3 
+ .06981141739 a; 2 - .001534877123 x + .009038433676 
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comparison of the experimental and approximated profiles 
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The first picture presents the experimental data, the second results from the 
fit-approximation (eqO) and the transition to the dimensionless frequencies . 

Now let's find the group-delay dispersion (GDD) in the active medium from 
the measured data: 

>#crystal length is 0.8 cm for double pass 

>GDD_cry_x := readdata('gdd_crystal_x.dat',l, float): 

>GDD_cry_y := readdata('gdd_crystal_y.dat', l,float): 

>n := vectdim(GDD_cry_y): 

>figl := plot([[GDD_cry_x[k],0.8*GDD_cry_y[k]] $k=l..n], 
color=red,axes=boxed) : 

>display(figl,title='measured GDD in active crystal'); 

>P1 := array(l..n): 
>Q1 := array(l..n): 

>for i from 1 to n do 

>Pl[i] := 
evalf((2*Pi*3*lelO/(GDD_cry_x[i]*le-7)-X_max)/bandwidth):# transition 
from wavelength to frequency 

>od: 

>fl := plot([[Pl[k],0.8*(bandwidthno- 15 )**2*GDD_cry_y[k]] 
$k= 1 . .n] ,color=red) : 

>eql := fit[leastsquare[[x,y], y=a*x 2 +b*x+c]]([[Pl[k] $k=l..n], 
evalf(0.8*(bandwidth*10 _15 )**2)*GDD_cry_y]);# fit-approximation 

>Q1 := [subs(x=Pl[k],rhs(eql)) $k=l..n]: 

>f2 := plot([[Pl[k],Ql[k]] $k=l..n],color=blue,axes=boxed): 

>display(fl,f2,title='GDD fit-approximation in frequency domain'); 
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measured GDP in active crystal 
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eql :=y= 5.557558914 x 2 + 21.42668994 x + 53.71773838 
GDP fit-approximation in frequency domain 
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eql is the result of the fit-approximation. 

The main devices for the GDD manipulation are the chirp mirrors. Their 
parameters are: 
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>#Chl (double pass) 

>GDD_Chl_x := readdata('gdd_chl_x.dat',l,float): 

>GDD_Chl_y := readdata('gdd_chl_y.dat',l,float): 

>n := vectdim(GDD_Chl_y): 

>plot([[GDD_Chl_x[k],2*GDD_Chl_y[k]] $k=l..n], 
color=red,axes=boxed,title='measured GDD for Chi'); 

>P2 := array(l..n-4): 
>Q2 := array (l..n-4): 

>for i from 1 to n-4 do 
>P2[i] := 
evalf((2*Pi*3*lelO/(GDD_Chl_x[i]*le-7)-X_max)/bandwidth):#transition 
from wavelength to frequency 

>Q2[i] := 2*(bandwidth*l<r 15 )**2*GDD_Chl_y[i] od: 

>fl := plot([[P2[k],Q2[k]] $k=l..n-4],color=red): 

>eq2 := fit[leastsquare[[x,y], 
y=a*x 6 +b*x 5 +c*x 4 +d*x 3 +e*x 2 +f*x+g]]([[P2[k] $k=l..n-4], [Q2[k] $k=l..n-4]]); 

>Q2 := [subs(x=P2[k],rhs(eq2)) $k=l..n-4]: 

>f2 := plot([[P2[k],Q2[k]] $k=l..n-4],color=blue): 

>display(fl,f2,axes=boxed,title='fit-approximation in frequency domain'); 
measured GDD for Ch1 
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eq2 := y = -4487.256771 x 6 + 8899.324767 x 5 - 4466.632623 x 4 - 612.3714657 x 3 
+ 692.5781343 x 2 - 17.99769657 x - 37.66103148 
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: 
-50 : 
-100 
-150 : 
-200 : 
-250 
-300 



fit-approximation in frequency domain 




-0.4 -0.2 0.2 0.4 0.6 0.! 



>#Ch2 (four passes) 

>GDD_Ch2_x := readdata('gdd_ch2_x.dat',l,fioat): 
>GDD_Ch2_y := readdata('gdd_ch2_y.dat',l,float): 
>n := vectdim(GDD_Ch2_y): 

>plot([[GDD_Ch2_x[k],4*GDD_Ch2_y[k]] $k=l..n], 
color=red,axes=boxed,title='measured GDD for Ch2'); 
>P3 := array(l..n): 
>Q3 := array(l..n): 

>for i from 1 to n do 
>P3[i] := 
evalf((2*Pi*3*lelO/(GDD_Ch2_x[i]*le-7)-X_max)/bandwidth):#transition 
from wavelength to frequency 
>od: 
>fl := plot([[P3[k],4*(bandwidth*l(T 15 )**2*GDD_Ch2_y[k]] 
$k= 1 . .n] ,color=red) : 

>eq3 := fit[leastsquare[[x,y], 
y=a*x 6 +b*x 5 +c*x 4 +d*x 3 +e*x 2 +f*x+g]]([[P3[k] $k=l..n], 
evalf(4*(bandwidth*lCT 15 )**2)*GDD_Ch2_y]); 
>Q3 := [subs(x=P3[k],rhs(eq3)) $k=l..n]: 
>f2 := plot([[P3[k],Q3[k]] $k=l..n],color=blue): 
>display(fl,f2,axes=boxed,title='fit-approximation in frequency domain'); 
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measured GDD for Ch2 
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eqS :=y= -2051.277573 x 6 + 4260.310370 x 5 - 3057.192323 x 4 + 565.9468050 x 3 
+ 385.5186662 x 2 - 74.89842760 x - 73.37050894 

fit-approximation in frequency domain 




And, at last, the GDD in output coupler and high-reflective mirrors: 

>#Oc (single pass) 

>GDD_Oc_x := readdata('gdd_oc_x.dat',l,float): 
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>GDD_Oc_y := readdata('gdd_oc_y.dat',l, float): 
>n := vectdim(GDD_Oc_y): 
>plot([[GDD_Oc_x[k],GDD_Oc_y[k]] $k=l..n], 
color=red,axes=boxed,title='measured GDD for OC); 
>P4 := array(L.n): 
>Q4 := array(l..n): 

>for i from 1 to n do 
>P4[i] := 
evalf((2*Pi*3*lelO/(GDD_Oc_x[i]*le-7)-X_max)/bandwidth):#transition 
from wavelength to frequency 
>od: 
>fl := plot([[P4[k],GDD_Oc_y[k]*(bandwidth*10- 15 )**2] $k=l..n], 
color=red): 

>eq4 := fit[leastsquare[[x,y], 
y=a*x 6 +b*x 5 +c*x 4 +d*x 3 +e*x 2 +f*x+g]] 

([[P4[k] $k=l..n], GDD_Oc_y*evalf((bandwidth*10- 15 )**2)]); 
>Q4 := [subs(x=P4[k],rhs(eq4)) $k=l..n]: 
>f2 := plot([[P4[k],Q4[k]] $k=l..n],color=blue): 
>display(fl,f2,axes=boxed,title='flt-approximation in frequency domain'); 

measured GDD for OC 
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eq4 :=y = 40.56991627 x 6 + 226.5828840x 5 - 15.88361104a; 4 - 6.140885675a; 3 
+ 1.492992122 x 2 + 1.204668103 x- .1174906766 
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fit-approximation in frequency domain 
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>#HR (four passes) 

>GDD_Hr_x := readdata('gdd_hr_x.dat',l, float): 
>GDD_Hr_y := readdata('gdd_hr_y.dat',l, float): 
>n := vectdim(GDD_Hr_y): 

>plot([[GDD_Hr_x[k],GDD_Hr_y[k]] 
$k=l..n],color=red,axes=boxed,title='measured GDD for HR'); 
>P5 := array(l..n): 
>Q5 := array(l..n): 
>for i from 1 to n do 
>P5[i] := 
evalf((2*Pi*3*lelO/(GDD_Hr_x[i]*le-7)-X_max)/bandwidth):#transition 
from wavelength to frequency 
>od: 
>fl := plot([[P5[k],GDD_Hr_y[k]*4*(bandwidth*l(T 15 )**2] 
$k= 1 . .n] ,color=red) : 

>eq5 := fit[leastsquare[[x,y], 
y=a*x 6 +b*x 5 +c*x 4 +d*x 3 +e*x 2 +f*x+g]]([[P5[k] 
$k=l..n], GDD_Hr_y*evalf(4*(bandwidth*l(n 15 )**2)]); 
>Q5 := [subs(x=P5[k],rhs(eq5)) $k=l..n]: 

>f2 := plot([[P5[k],Q5[k]] $k=l..n],color=blue): 
>display(fl,f2,axes=boxed,title='fit-approximation in frequency domain'); 
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measured GDD for HR 
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eq5 :=y = 371.6378105 a; 6 + 1859.815396 x 5 - 146.0722578 x 4 - 186.9742695 a; 3 
+ 12.76215250 x 2 + 13.36697439 x- .5561160074 

fit-approximation in frequency domain 




Now, as a result of the obtained fit-approximations, we have the normalized 
net-GDD with corresponding FORTRAN-code for simulation: 

>end res := 
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evalf(rhs(eql)+rhs(eq2)+rhs(eq3)+rhs(eq4)+rhs(eq5)); 

>plot(%,x=-0.45..0.8,axes=boxed,title='normalized to tf net-GDD'); 
>codegen[fortran] (end_res) ; 

end_res := -57.98740873 - 6126.326617 x 6 + 15246.03342 x 5 - 7685.780815 x 4 
- 239.5398159 x 3 + 1097.909504 x 2 - 56.89779174 x 

normalized to tf net-GDD 
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tO = 
-0 . 5798741E2-0 . 6126327E4*x**6+0 . 1524603E5*x**5-0 . 7685781E4*x* 

#*4-0 . 2395398E3*x**3+0 . 109791E4*x**2-0 . 5689779E2*x 
Or in the usual co-ordinates: 

>P6 := array(1..100): 

>Q6 := array(1..100): 

>for i from 1 to 100 do 

>P6[i] := 2*Pi*3elO*le7/(bandwidth*(-0.4+0.8*i/100)+X_max):# from 
frequency to wavelength [in nm] 

>Q6[i] := 
evalf(subs(x=-0.4+0.8*i/100,end_res/(bandwidth*le-15) 2 )): 

>od: 

>plot([[P6[k],Q6[k]] $k=1..100], 
color=red,axes=BOXED,title='net-GDD [fs 2 ] vs. wavelength [nm]'); 
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net-GDD [fs A 2] vs. wavelength [nm] 
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The experimentally observed ultrashort-pulse spectrum is: 

>Spec_x := readdata('LiSGaF_x.dat',l,float): 
>Spec_y := readdata('LiSGaF_y.dat',l, float): 
>n := vectdim(Spec_y): 

>plot([[Spec_x[k],Spec_y[k]] 
$k=l..n],color=red,axes=boxed,title='experimental spectrum (a.u. vs. wave- 
length)'); 

>P5 := array(l..n): 
>Q5 := array(l..n): 
>for i from 1 to n do 
>P5[i] := 
evalf((2*Pi*3*lelO/(Spec_x[i]*le-7)-X_max)/bandwidth):#transition 
from wavelength to frequency 
>Q5[i] := Spec_y[i]: 
>od: 

>plot([[P5[k],Q5[k]] $k=l..n], 
color=green,axes=boxed,title='experimental spectrum (a.u. vs/ dimensionless 
frequency)'); 
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experimental spectrum (a.u. vs. wavelength) 
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experimental spectrum (a.u. vs/ dimensionless frequency) 
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We can see a pronounced red shift relatively to the gain band 
center. The narrow line is the so-called side-band or dispersive 
wave. 



The typical results of the numerical simulations based on the described above 
model are presented in the next figure (power is given in the arbitrary units). 

>Specl_x := readdata('numl_x.dat',l, float): 
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>Specl_y := readdata('numl_y.dat',l,float): 
>Spec2_x := readdata('num2_x.dat',l, float): 
>Spec2_y := readdata('num2_y.dat',l, float): 

>Spec3_x := readdata('num3_x.dat',l, float): 
>Spec3_y := readdata('num3_y.dat',l, float): 
>Spec4_x := readdata('num4_x.dat',l, float): 
>Spec4_y := readdata('num4_y.dat',l,float): 
>Spec5_x := readdata('num5_x.dat',l, float): 
>Spec5_y := readdata('num5_y.dat',l, float): 

>Spec6_x := readdata('num6_x.dat',l, float): 
>Spec6_y := readdata('num6_y.dat',l, float): 
>n := vectdim(Specl_y); 
>pl := plot([[Specl_x[k],Specl_y[k]] $k=l..n/4],color=black): 
>n := vectdim(Spec2_y); 

>p2 := plot([[Spec2_x[k],Spec2_y[k]] $k=3*(n-l)/4..n], 
color=black): 

>n := vectdim(Spec3_y); 
>p3 := plot([[Spec3_x[k],2.5*Spec3_y[k]] $k=3*(n-l)/4..n], 
color=red): 

>n := vectdim(Spec4_y); 
>p4 := plot([[Spec4_x[k],2.5*Spec4_y[k]] $k=l..n/4], 
color=red): 

>n := vectdim(Spec5_y); 

>p5 := plot([[Spec5_x[k],25*Spec5_y[k]] $k=3*(n-l)/4..n], 
color=blue): 

>n := vectdim(Spec6_y); 

>p6 := plot([[Spec6_x[k],25*Spec6_y[k]] $k=l..n/4], 
color=blue): 

>display(pl,p2,p3,p4,p5,p6,view=0..1,axes=BOXED,title='numerical spec- 
tra (a.u.)'); 
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numerical spectra (a.u. 
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In the "natural" co-ordinates we have: 

>n := vectdim(Specl_y); 

>pl:=plot([[2*Pi*3elO*le7/(bandwidth*Specl_x[k]+X_max),Specl_y[k]] 
$k=l..n/4],color=black): 

>n := vectdim(Spec2_y); 
>p2:=plot([[2*Pi*3el0*le7/(bandwidth*Spec2_x[k]+X_max),Spec2_y[k]] 
$k=3*(n-l)/4..n],color=black): 
>n := vectdim(Spec3_y); 
>p3 := 
plot([[2*Pi*3el0*le7/(bandwidth*Spec3_x[k]+X_max),2.5*Spec3_y[k]] 
$k=3*(n-l)/4..n],color=red): 

>n := vectdim(Spec4_y); 
>p4 := 
plot([[2*Pi*3el0*le7/(bandwidth*Spec4_x[k]+X_max),2.5*Spec4_y[k]] 
$k=l..n/4],color=red): 

>n := vectdim(Spec5_y); 
>p5 := 
plot([[2*Pi*3el0*le7/(bandwidth*Spec5_x[k]+X_max),25*Spec5_y[k]] 
$k=3* (n-1) /4. .n] ,color=blue) : 
>n := vectdim(Spec6_y); 
>p6 := 
plot([[2*Pi*3el0*le7/(bandwidth*Spec6_x[k]+X_max),25*Spec6_y[k]] 
$k=l..n/4],color=blue): 

>display(pl,p2,p3,p4,p5,p6,view=0..1,axes=BOXED,title='pulse spectrum 
vs. wavelength [nm]'); 
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julse spectrum vs. wavelength [nm] 




750 



800 



850 



900 



950 



1000 



The black curve corresponds to the soliton propagation without lasing fac- 
tors, i.e. in the absence of gain, loss saturation and spectral filtering. The 
high-order dispersion slightly transforms the spectrum of the soliton, but there 
is not the visible frequency shift. Perhaps, the soliton duration (169 fs) is too 
large for the high-order dispersions' manifestation. But we can not change the 
pulse duration for fixed GDD. Such possibility is opened by lasing in the pres- 
ence of the gain and saturable loss. The red and blue curves correspond to the 
above described laser model with 1.2 W absorbed pump power for 20x30 /i 2 
pump mode and gaussian mode with 25 fim diameter, correspondingly. The 
pulse durations in these cases are 31 and 27 fs, correspondingly. We can see the 
appearance of spectral spikes (dispersion waves), which locate near from zeros 
and 2 7rm - values of GDD. The domain of the large gradient of GDD in the red 
spectral region is filled by such spikes. But the pulses' spectra have the very 
small frequency shifts from gain band maximum, although there is the spectral 
"shoulder" in the red region. The last corresponds to local GDD extremum near 
from 883 nm (see above). 



We had performed the various manipulations with the pump 
and the fast absorber saturation intensity, but as a result, the 
self-frequency shift can not be obtained only due to high-order 
dispersion. 



So, there is some additional mechanism of the pulse spectrum shift. As 
the pulse duration is greater than 10 fs, the nonlinear dispersion can not cause 
such shift. Therefore we have to investigate the stimulated Raman scattering 
influence on the pulse spectrum. 
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Let assume, that the frequency shift results from the stimulated Raman 
scattering in the active crystal. The vibrational amplitude Q in the dependence 
on the pump and signal spectral amplitudes A_ m and A_ n, respectively, obeys 
the following equation (see |l6|): 

>restart: 

>eql := 
diff(Q(t),t$2)+2*diff(Q(t),t)/T+Omega 2 *Q(t) = 
mu*A_m*exp(I*m*omega*t)*conjugate(A_n)*exp(-I*n*omega*t); 

>Q = subs({_Cl=0,_C2=0},rhs(dsolve(eql,Q(t)))); 

>amp_Q = 
Sum(Sum(expand(numer(rhs(%))/expand(denom 
(rhs(%)))/exp(-I*t*omega*(-m+n))),n=-N/2..N/2),m=-N/2..N/2); 
#here we supposed that the pulse width « T (T is the phonon relaxation time) 



eql :=(^Q(t))+ 2{ ^ {t)) +n 2 Q(t)= f ,A_ m ^ Im -^TAjr)e^ In -^ 



liA_me { ~ IuJt{m - n)) {A_n)T 



Tfl 2 + 2Imuj — 2 I nuj — T lo 2 m 2 + 2 Tlu 2 mn — T lo 2 n 2 



amp _ 



nA_m(A_n)T 



1/2N I 1/2 iV 

^ 1 <H T n 2 + 2 1 ma; - 2 1 noj - T u 2 m 2 + 2T u 2 mn - T cu 2 n 2 

=-1/2 JV \n=-l/2JV 



Here m andn are the mode numbers, w is the frequency interval between 
field spectral components, fl is the Raman frequency, T is the relaxation time, 
N is the number of the frequency components in the field spectrum (in the case 
of our numerical simulations N= 2 13 ), jj, is the positive real number expressing 
photon-phonon coupling. Note, that the last expression for the amplitude of 
the vibrational oscillations (the frequency of these oscillations is (m-n)* u>) can 
be re-written as 

>amp_Q = 

Sum(Sum(mu*A_m*conjugate(A_n)/ (Omega 2 - 
2*I*omega*(m-n)/Tr - omega 2 *(n-m ) 2 ), 
n = -N/2 .. N/2),m = -N/2 .. N/2); 

>amp_Q = 

Sum(Sum(mu*A_m*conjugate(A_n)*(Omega 2 - omega 2 *(n-m) 2 + 
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2*I*omega*(m-n)/Tr)/((Omega 2 - omega 2 * (n-m) 2 ) 2 + 
4*omega 2 *(m-n) 2 /(Tr 2 )),n = -N/2 .. N/2),m = -N/2 .. N/2); 



1/2 JV / 1/2 JV -r-j r 

jjl A_m \A_n) 



amp_Q= ^2 ^2 

m=-l/2JV \n=-l/2JV fj 2 -| — ~^"" — - u 2 (-171 + n) 2 I 

\ TV ' 



m=-l/2N 



2Iu) (to — n) 

l 2N /j,A_m(A_n)(n 2 -oj 2 (-m + n) 2 + 



i'-v ( 1/2 JV Uj4 m U nU» 2 -^ 2 f-TO + n) 2 + 2/U;(m n) ^ 

Tr 



amp V = / 7 ,.>■■> 



. n=-l/2 JV (fi2 _ ^2 (_ m _|_ n )2)2 _| _ 

\ TV 



Now we have to define the values of the used parameters. There are three 
Raman lines in LiSGaF with following characteristics: fi=551, 349, 230 cm^ 1 ^, 
1/T = 6.2, 7.6, 4.2 cm^ 

>Omegal := evalf(2*Pi*3*10 1 0*551):#[Hz] 
>Omega2 := evalf(2*Pi*3*10 1 0*349):#[Hz] 
>Omega3 := evalf(2*Pi*3*10 1 0*230):#[Hz] 

>bandwidth := .481859640el5:#[Hz] the normalization for frequen- 
cies 

>Omegal := evalf(Omegal/bandwidth);#normalized Omegal 
>Omega2 := evalf(Omega2/bandwidth);#normalized Omega2 
>Omega3 := evalf(Omega3/bandwidth);#normalized Omega3 
>omega=evalf(2*Pi*bandwidth/N):#[Hz] 

>print('normalized w:'); 
>2*Pi/N;#normalyzed omega 
>T1 := evalf(l/(2*Pi*3*10 1 0*6.2)):#[s] 
>T2 := evalf(l/(2*Pi*3*10 1 0*7.6)):#[s] 
>T3 := evalf(l/(2*Pi*3*10 1 0*4.2)):#[s] 

>T1 := evalf(Tl*bandwidth);#normalized relaxation time 

>T2 := evalf(T2*bandwidth);#normalized relaxation time 
>T3 := evalf(T3*bandwidth);#normalized relaxation time 
>solve(gain_s=6*omega_s*chi/n_s/c,chi):#if permittivity chi in [cm 2 /W], 
gain_s is the Raman signal gain 

>mu := simplify(3*omega_s*%*2*Omega/n_s/c/T); 

>mul := subs({Omega = Omegal,T = Tl,gain_s=evalf(2.5*1.2e-10), 
beta=.3379192098e-ll},mu*0.8/beta);#normalized 

>mu2 := subs({Omega = Omega2,T = T2,gain_s=evalf(. 185*1. 2e-10), 
beta=.3379192098e-ll},mu*0.8/beta);#normalized 

>mu3 := subs({Omega = Omega3,T = T3,gain_s=evalf(.15*1.2e-10), 
beta=.3379192098e-ll},mu*0.8/beta);#normalized 

m := .2155421299 
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02 := .1365230551 



03 := .08997221396 



normalized w 



2^ 
TV 



Tl := 412.3136751 



T2 := 336.3611560 



T3 := 608.6535205 



gain s O 
V '■= m 



/il := .03712810569 



[il := .002133193454 



H'S := .0006299236009 

So, we can investigate the field evolution on the basis of the following equa- 
tion: 

>Diff(A_n(z),z) = 
Sum(mu*conjugate(A_m)*A_n*(I*(Omega 2 - omega 2 *(n - m) 2 )+2*omega*(m 
- n)/Tr)/( (Omega 2 - omega 2 *(n - m) 2 ) 2 + 4*omega 2 *(m - n) 2 /(Tr 2 ))*A_m,m 
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= -N/2 .. N/2); 

>Diff(A_n(z),z) = 
A_n*Sum(mu*abs(A_m) 2 *(I*(Omega 2 - omega 2 *(n - m) 2 ) + 2*omega*(m - 
n)/Tr)/( (Omega 2 - omega 2 *(n - m) 2 ) 2 + 4*omega 2 *(m - n) 2 /(Tr 2 )),m = -N/2 
•• N/2); 



d_ 
dz 



A_n(z) 



V2W gam 8 Sl(A m)A n(I (ft 2 - uj 2 (-m + n) 2 ) + 1^S™ — Vl) A m 

E " " [r 

m=-l/2iV 



T ((ft 2 - uo 2 (-to + n) 2 ) 2 + 



.,.., 4lu 2 (to - n) 2 



Tr z 



— A_n(z) = 
oz 



A_n 


( i/2W gam sftM ml 2 (/(ft 2 cv 2 ( m + n) 2 ) + 2u ' ^ 
ST^ ' ' Tr 


n)A 




^^ //„o n , n^no 4o> 2 (m — n) 2 . 

.m— 1/2 JV T((ft 2 w 2 ( TO + n) 2 ) 2 + 2 ) 

\ Tr 


/ 



Here n is the "signal" mode, m is the "pump" mode, degenerate case m=n 
corresponds to the pure self-phase modulation, Stokes signal shift m>n cor- 
responds to an amplification, anti-Stokes shift m<n corresponds to a loss of 
signal wave. Because of the Raman line is narrow in the comparison with pulse 
spectrum (see normalized T), we can re- write the expressions: 

>Diff(A_n(z),z) = 
I*A_n*Sum(mu*abs(A_m) 2 /(Omega 2 + 2*I*omega*(m - n)/Tr-omega 2 *(n - 
m) 2 ),m = -N/2 .. N/2); 

>Diff(A_n(z),z) = 
I*A_n*Sum(mu*abs(A_m) 2 /((Omega - omega*(n - m))*(Omega + omega*(n 
-m )) + 2*I*omega*(m - n)/Tr),m = -N/2 .. N/2); 

>Diff(A_n(z),z) = 
I*A_n*Sum(mu*abs(A_m) 2 /(2*Omega*(Omega - omega*(n - m)) + 
2*I*0mega/Tr),m = -N/2 .. N/2); 



/ i 



-—A n(z) = I A n 
oz 



/2N 

E 



gain _s ft |^4_ni|' ! 



™ ,ni o/ xo 2Iu(m — n), 

Vm=-l/2 N T (ft 2 - U 2 (-TO + n) 2 H 7- 

\ Tr 
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9 A , x 

[ v-^ nam s 11 \A m\ 

IA_n \ Yl ■ - I - 



„//„ , xx /„ / xx 2 1 u (m — n) N 
m=-i/2JvT((fi - w(-m + n)) (fi + w(-m + n)) H - -), 



■I I 1 / 2N ■ O I A I 2 

-A_n(*) = .M_n £ " " 



o JA - xx V^ / ^^ A ' * I / ^ 9 A O 

\m=-i/2;vT(2fi(fi-u;(-m + n)) + V 

\ Tr 



The last expression results from the assumption T»l. In this case 

>Diff(A_n(z),z) = 
I*A_n*mu*abs(A_m) 2 *Int((2*Omega*(Omega - x) - 2*I*0mega/Tr)/ 
(4*Omega 2 * (Omega - x) 2 + 4*Omega 2 /(Tr 2 )),x = -infinity .. infinity); 



oo 2 \i 


(0- 


x) + 


-2in 

Tr 


dx 


-°°4tt 


2 (Q- 


-x) 2 


4S1 2 

H 9 



I A_ngain_sfl \A_m\ 

— A n(z) = 

dz ~ y ' T 



In the last expression x is the frequency difference, A_m corresponds to 
the pump intensity for - w(n-m)= ft. The contribution of the real part in the 
integral (the self-phase modulation from the both sides of the Raman line) is 
equal to 0. 

>-(I*Tr/2/Omega)*Int(l/(x 2 +l),x=0.. infinity); 
> value (%); 



-1 f°° 1 

— I Tr — dx 

2 Jo x 2 + l 



I Trir 

_4 

n 

Thus we obtained the simplest expressions for the field evolution: 
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>Diff(A_n(z),z) = A_n*gain*Pi*abs(A_m) 2 /4; 
>Diff(A_m(z),z) = -A_m*gain*Pi*abs(A_n) 2 /4; 



d 1 2 

——A n(z) = —A n gaimr \A ml 
oz 4 



d 1 ,2 

— A m(z) = — A m gainn \A n\ 
oz 4 



The stimulated Raman gain parameters for numerical simulation are: 

>gain_sl := evalf(0.8*2.5*1.2e-10/.3379192098e-ll); 
>gain_s2 := evalf(0.8*. 185*1. 2e-10/.3379192098e-ll); 
>gain_s3 := evalf(0.8*.15*1.2e-10/.3379192098e-ll); 



gain_sl := 71.02289335 



gain_s2 := 5.255694108 



gain_s3 := 4.261373601 



Additionally we are to take into consideration the spontaneous Raman scat- 
tering as a source for the stimulated scattering. The gain coefficients is this case 
are @: 

>d_sigma := 
(4*omega_l/3/omega_s)*3*n_s 2 *h*omega_s 3 * 
gain_s/(8*Pi 3 *c 2 *N)/(l - exp(-h*(omega_l - 

omega_s)/(2*Pi*kb*Tc)));#kb and Tc are the Boltzmann's constant and tem- 
perature, correspondingly 

>d_sigmal := 
evalf(sqrt(subs( 

{kb=1.38*le-23,Tc=300,omega_s=2*Pi*3el0/.9051e-4, 
omega_l=2*Pi*3el0/.85e-4,n_s=1.4,h=6.62e-34, 
gain_s=gain_sl,c=3elO,N=le20},d_sigma*N))); 

>d_sigma2 := 
evalf(sqrt(subs( 
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{kb=1.38*le-23,Tc=300,omega_s=2*Pi*3el0/.88e-4, 

omega_l=2*Pi*3el0/.85e-4,n_s=1.4,h=6.62e-34, 

gain_s=gain_s2,c=3elO,N=le20},d_sigma*N))); 

>d_sigma3 := 
evalf(sqrt(subs( 

{kb=1.38*le-23,Tc=300,omega_s=2*Pi*3el0/.87e-4, 
omega_l=2*Pi*3el0/.85e-4,n_s=1.4,h=6.62e-34, 
gain_s=gain_s3,c=3elO,N=le20},d_sigma*N))); 

1 omega _l omega _s 2 n_s 2 h gain _s 
a siqma : = ■ —, i ; — 

_ il 9 O T »r /, I 1 /n h (omega _l-omega_s) . 

^ tt a c z N (I — e ( l/A T,kbTc >) 

d_sigmal := .0001280998590 

d_sigma2 := .00003815471494 

d_sigma3 := .00003767032764 



d_ sigma are the increments of the spontaneous Stokes components (i.e. Raman 
spontaneous seeds) growth . 

As a result of the simulations on the basis of this model, we obtained next 
spectra: 

>with('linalg'): 
>with(plots): 

>R_Specl_x := readdata('specl_x.dat',l, float): 
>R_Specl_y := readdata('specl_y.dat',l, float): 

>R_Spec2_x := readdata('spec2_x.dat',l,float): 
>R_Spec2_y := readdata('spec2_y.dat',l, float): 
>R_Spec3_x := readdata('spec3_x.dat',l, float): 
>R_Spec3_y := readdata('spec3_y.dat',l,float): 
>R_Spec4_x := readdata('spec4_x.dat',l, float): 
>R_Spec4_y := readdata('spec4_y.dat',l, float): 

>R_Spec5_x := readdata('spec5_x.dat',l,float): 
>R_Spec5_y := readdata('spec5_y.dat',l, float): 
>R_Spec6_x := readdata('spec6_x.dat',l, float): 
>R_Spec6_y := readdata('spec6_y.dat',l,float): 
>n := vectdim(R_Specl_y); 

>pl := plot([[R_Specl_x[k],R_Specl_y[k]] 
$k= 1 . .n/ 16] ,color=black) : 
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>n := vectdim(R_Spec2_y); 
>p2 := plot([[R_Spec2_x[k],R_Spec2_y[k]] 
$k=7*(n-l)/8..n],color=black): 

>n := vectdim(R_Spec3_y); 
>p3 := plot([[R_Spec3_x[k],R_Spec3_y[k]] 
$k= 1 . .n/ 16] ,color=red) : 

>n := vectdim(R_Spec4_y); 

>p4 := plot([[R_Spec4_x[k],R_Spec4_y[k]] 
$k=7*(n-l)/8..n],color=red): 

>n := vectdim(R_Spec5_y); 
>p5 := plot([[R_Spec5_x[2*k],R_Spec5_y[2*k]] 
$k=l..n/32],color=blue): 

>n := vectdim(R_Spec6_y); 
>p6 := plot([[R_Spec6_x[2*k],R_Spec6_y[2*k]] 
$k=7*(n-l)/16..n/2],color=blue): 

>display(pl,p2,p3,p4,p5,p6,axes=boxed, title='pulse spectrum vs. dimen- 
sionless frequency'); 

pulse spectrum vs. dimensionless frequency 




-0.4 -0.3 -0.2 -0.1 0.1 0.2 



>X_max := .2245909122el6: 
>bandwidth := .481859640el5: 
>n := vectdim(R_Specl_y): 
>pl := 
plot([[2*Pi*3elO*le7/(bandwidth*R_Specl_x[k] + X_max),R_Specl_y[k]] 
$k= 1 . .n/ 16] ,color=black) : 

>n := vectdim(R_Spec2_y): 
>p2 := 
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plot([[2*Pi*3elO*le7/(bandwidth*R_Spec2_x[k] + X_max),R_Spec2_y[k]] 
$k=7*(n-l)/8..n],color=black): 

>n := vectdim(R_Spec3_y): 
>p3 := 
plot([[2*Pi*3el0*le7/(bandwidth*R_Spec3_x[k] + X_max),R_Spec3_y[k]] 
$k= 1 . .n/ 16] ,color=red) : 

>n := vectdim(R_Spec4_y): 
>p4 := 
plot([[2*Pi*3el0*le7/(bandwidth*R_Spec4_x[k] + X_max),R_Spec4_y[k]] 
$k=7*(n-l)/8..n],color=red): 

>n := vectdim(R_Spec5_y): 
>p5 := 
plot([[2*Pi*3el0*le7/(bandwidth*R_Spec5_x[2*k] + X_max),R_Spec5_y[2*k]] 
$k=l..n/32],color=blue): 

>n := vectdim(R_Spec6_y): 
>p6 := 
plot([[2*Pi*3el0*le7/(bandwidth*R_Spec6_x[2*k] + X_max),R_Spec6_y[2*k]] 
$k=7*(n-l)/16..n/2],color=blue): 

>display(pl,p2,p3,p4,p5,p6,axes=BOXED, 
title='pulse spectrum vs. wavelength [nm]'); 

pulse spectrum vs. wavelength [nm] 
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Three different spectra are presented in the figure. Black curve corresponds 
to the soliton propagation (200 fs pulse duration) in the presence of stimulated 
Raman scattering. As a result of the propagation, the Stokes component appears 
at the main Raman frequency. The strong scattering destroys the soliton after 10 
000 cavity transitions. In the laser, a balance between all lasing factors stabilizes 
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the ultrashort pulse. But there is a visible red shift of the field spectrum (red 
and blue curves) due to stimulated Raman scattering. The pulse spectrum can 
be pushed from the gain band center (red curve, 25 fs pulse duration) as a result 
of the Raman self-scattering. There is the possibility of the generation of the 
additional Stokes lines (blue curve, 56 fs). 



The frequency shift is comparable with the experimental one. 
Hence, the stimulated Raman scattering is the main source of 
the Stokes shift of the pulse spectrum in the Cr:LiSGaF Kerr- 
lens mode-locked laser (see |lq|). 



The collaboration between Maple and external numerical simulators (based 
on the FORTRAN-code in our case) proves to be extremely fruitful: 



1) analytical model building (Maple) =4> 2) external code gen- 
eration (Maple) =>■ 3) calculation of the simulation parameters 
(Maple) =>■ 4) external calculations (FORTRAN-code, the exter- 
nal program can be started from the Maple directly through the 
"system" call.) =>■ 5) data processing (Maple) =4> 6) analytical 
interpretation of the results (Maple) . 



As an example of the last step, let's consider the problem of the pulse sta- 
bility in the Kerr-lens mode-locked laser. 



9.2 Multipulsing and ultrashort pulse stability 

We try to shorten the pulse duration and increase its energy. For this aim we 
tend the net-GDD to zero (see Part VIII). As a result, the ultrashort pulse 
stability can be lost. Let's consider this phenomenon in detail (||l£]). In the 
framework of the abberationless approximation the stability loss is revealed 
as the absence of the soliton-like as well as breezer solution (the evolutional 
equations for the pulse parameters diverge) . What is meaning of this divergence? 
The answer comes from the numerical simulations based on the above de- 
scribed model. Let's neglect the stimulated Raman scattering and the higher- 
order dispersions. In this case, the typical results of the simulations demonstrate 
multipulsing in the vicinity of zero GDD. 



Quasi-soliton consideration fails due to the appearance of the 
regular or irregular multiple pulse generation. 



The boundary of the soliton-like pulse nonstability obtained from the numerical 
simulation is shown here (the parameters in question can be found in the work 
of reference) : 

>restart: 
>with(plots): 
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>with(stats): 

>points_numer_x := 
[-160, -75, -40, -25, -18, -16, -20, -40, -45, -75, -150]:# GDD is normal- 
ized to tf 

>points_numer_y := [0.2, 0.5, 1, 2, 5, 10, 20, 30, 80, 200, 500] :# sigma is nor- 
malized to the self-phase modulation coefficient beta 

>statplots[scatterplot](points_numer_x, evalf( 
map(logl0,points_numer_y) ),axes=boxed,color=red,symbol=box): 

>display(%,color=red,TEXT([-120,l], "stable single pulse")): 

>figl := 
display ( 

%,color=blue,TEXT([-40,2.6],"multipulsing"),TEXT([-50,-0.6], 
"multipulsing"),title='logarithm of boundary sigma vs. GDD'): 



> display (figl); 



logarithm of boundary sigma vs. GDD 



1-U^- 



multipulsing 



stable single pulse 



multipulsing 



-160 -140 -120 -100 -80 -60 -40 -20 



There are the upper and lower on a boundaries defining the 
transition to multipulsing. There is no single pulse in the vicin- 
ity of zero GDD. 



To interpret this picture let's return to the generalized Landau- Ginzburg 
equation in the week-nonlinear limit, when o $> <Cl (see Part VII). For the 
steady-state ultrashort pulse propagation we have (time and intensity $ are 
normalized): 
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>master_l := = alpha - gamma + I*phi + diff(rho(t),T(t,2))/rho(t) 
+ I*k_2*diff(rho(t),T(t,2))/rho(t) + gamma*sigma*Phi(t) - I*Phi(t);# rho is 
the filed, Phi is the intensity 

>fl := (t)->rhoO*sech(t*tau) 1+/ * psi ;# quasi-soliton profile 

>f2 := (t)->rho0 2 *sech(t*tau) 2 ;# pulse intensity 



jLp(t) Ik 2(jLp(t)) 
master _1 := = a - 7 + I (f>+ w - y ->- ~ y ^ yy — 



Pit) 



Pit) 



■ya$(t) - I$(t) 



fl :=t^ pOscch(tT) (1+I ^ 



f2 :=«-► p0 2 sech(tr) 



> simplify ( 
subs({rho(t)=fl(t),Phi(t)=f2(t)},rhs(master_l)) ): 
>numer(%): 



>eql 
>eq2 
>eq3 



collect(%,cosh(t*tau) 2 ): 

evalc( coeff(eql,cosh(t*tau),2) ); 

evalc( coeff(eql,cosh(t*tau),0) ); 



eq2 := a - 7 + t 2 - t 2 ip 2 - 2 k_2 t 2 ip + I {(j> + 2t 2 ip + k _2 t 2 - k_2 t 2 ip 2 ) 



2 „/,2 „ n 2 



eq3 :=japO z + T 2 ip z -2T z + 3 k_2 t z ip + I (-3 t z ip + k_2 t z ip z - pO - 2 k_2 t z ) 



>eq4 := coeff(eq2,I,0); 



eq4 := a - 7 + t z - t z ijr - 2 k_2 t z ip 



As it was shown in Parts VII, VIII, the pulse exists if a — 7 < 



Simul- 
taneously, it is condition of the cw suppression (the threshold is not exceeded 
for the noise out of pulse). As a result, the stability against cw oscillation is 
provided with: 



>eq5 := factor( eq4 - (alpha-gamma) )/tau 2 > 0; 
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eq5 :=0 < 1 - tp 2 -2k_2tp 



The pulse chirp can be found from the equation eq3: 

>eq6 := subs({rho0 2 =x,tau 2 =y},coeff(eq3,I)=0): 
>eq7 := subs({rho0 2 =x,tau 2 =y},coeff(eq3,I,0)=0): 

>simplify (subs( x=solve(eq7,x),eq6 ) ):#we find the intensity 
>numer(lhs(%))/y = 0: 
>sol := solve(%,psi); 

I 3 7 cr-3fc_^ + y / 97 2 cr 2 -2k_2ja + 9k_2 2 + 8 + 8k_2 2 j 2 a 2 
S ° l := 2 l+k_2 1 a ' 

1 3ja-3k_2- y / 9 7 2 cr 2 - 2k_2 ^a + 9 k_2 2 + 8 + 8k_2 2 ~/ 2 a 2 

2 l + k_2-/a 

In the combination with the condition eq5 we have: 

>eq8 := solve( numer( simplify ( subs(psi=sol[l],rhs(eq5)) ) ) = 
0, sigma ): 

>eq9 := solve( numer( simplify ( subs(psi=sol[2],rhs(eq5)) ) ) = 
0, sigma ): 

>plot([logl0(subs(gamma=0.01,eq8[l])),logl0(subs(gamma=0.01,eq8[2])), 
Iogl0(subs(gamma=0.01,eq9[l])),logl0(subs(gamma=0.01,eq9[2]))], 
k_2=-160..0,axes=boxed,color=[green, magenta]): 

>display(%,figl,view=[-160..0,-1..3]); 
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logarithm of boundary sigma vs. GDP 




Hence the lower boundary of the pulse destabilization is good 
approximated by the condition of the cw excitation (magenta 
curve) for the large negative GDD. The green curve corresponds 
to the chirp-free generation in the soliton model. The intersec- 
tion of this curve with the stability boundary (red points) gives 
the system's parameters corresponding to the minimal pulse du- 
ration. 



However, there are the differences between the analytical model and the numer- 
ical results: 1) the ultrashort pulse is chirp-free in the wider region than that 
predicted from the soliton model; 2) the a increase providing the pulse shorten- 
ing causes the pulse destabilization (upper on the a parameter boundary of the 
pulse stability, red points); 3) the pulse is unstable in the vicinity of zero GDD. 



The numerical simulations demonstrate, that the pulse destabi- 
lization on the upper stability boundary occurs for the negative 
net-gain coefficient a — 7 <0. This prevents the cw excitation. 



Let's consider the toy model of the pulse destabilization in the absence of 
the cw excitation. For the sake of simplicity, we shall consider the totally real 
Landau- Ginzburg equation (Part VI). 

>master_2 := rho(t)*g + diff(rho(t),T(t,2)) + 
rho(t) 3 *Sigma;#g=alpha-gamma, Sigma=gamma*sigma 
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master _2 := p(t) g + (— p(i)) + p(£) 3 £ 



Let's expand this equation on the small perturbation £(£): 

>expand( subs( rho(t)=rho(t)+mu*zeta(t),master_2 )): 
>limit( diff(%,mu),mu=0 ):# functional Frechet derivative 
>master 3 := master 2 + %; 



d 2 d 2 

master _3 := p{t) g + (^ /»(*)) + pW' 2 + .9 C(*) + (g^ C(*)) + 3 S p(t) 2 <(i) 



and find its steady-state solutions. Thereto we make the following substitution: 

>f3 := (t)->rhoO*sech(t*tau); 

>f4 := (t)->epsilon*diff( sech(t*tau) ,t$2);#rho0 



f3 :=£-*■ pOsech(ir) 



# :=i^e(^sech(tr)) 



> simplify ( 
subs({rho(t)=f3(t),zeta(t)=f4(t)},master_3) ): 
>expand( numer(%)/rhoO): 
>eqlO := collect ( numer(%),cosh(t*tau) ); 

eqlO := (pO 5 + g £ r 2 + pO t 2 + e r 4 ) cosh(t r) 4 

+ (-2 5 e r 2 + pO 3 S + 3 S pO 2 e r 2 - 2 pO r 2 - 20 £ t 4 ) cosh(t r) 2 + 24 e r 4 

- 6 E pO 2 e t 2 

We can see, that this substitution obeys the perturbed steady-state equation 
( £ is the perturbation amplitude). 

>eqll := expand( 
coeff(eqlO,cosh(t*tau),0)/epsilon/tau 2 ); 
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>eql2 := expand( coeff(eql0,cosh(t*tau),2) ); 
>eql3 := simplify ( coeff(eqlO,cosh(t*tau),4)); 



eqll :=24t 2 - 6 E pO 2 



eql2 := -2ger 2 + pO 3 E + 3 E pO 2 £ r 2 - 2 pOr 2 -20et 



eg!5 := pO .g + g e r 2 + pO t 2 + e r 4 



>sol := solve({eqll=0,eq!2=0,eql3=0},{rho0,tau,epsilon}); 



sol := {t = RootOi(g + _Z 2 ), pO = 2RootOf(E _Z 2 + g, label = _L1), 

2 RootOffE Z 2 +g, label = LI) , , 

e = = ■ = }, {e = e, pO = 0, t = 0}, { 

3 9 
2 

pO = - g RootOf (_Z 2 g E + 5, label = _L2), 
5 

e = 2 RootOf (_Z 2 g E + 5, label = _L2), 
t = RootOf (.g + 5 _Z 2 )} 

>soll_tau := allvalues(subs(sol[l],tau)); 
>soll_rho := allvalues(subs(sol[l],rhoO)); 

>soll_e := allvalues(subs(sol[l],epsilon)); 
>sol2_tau := allvalues(subs(sol[3],tau)); 

>sol2_rho := allvalues(subs(sol[3],rho0)); 
>so!2_e := allvalues(subs(sol[3],epsilon)); 



soil _tau := \f—g, —\f—g 



soil rho :=2,/-l -2 . ' 9 



9_ _9_ 

2 V E 2 V E 
soil _e := , 

3 9 3 g 
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sol2 tau := \ — 



1 



2/12/1 
sol2_rho:=- 9] /-5— ,-- 9] j-5 — 



soli < ■= 2 ( /-5— -, -2 a/-5— - 
gS V 5 s 



There exist two types of the perturbed solutions . The first one correspond- 



ing to unperturbed solution for the arbitrary small e has a form 



>soll := subs( 
{rhoO=soll_rho[l],tau=soll_tau[l],epsilon=soll_e[l]}, 
expand((f3(t) + f4(t)) 2 ) ): 

>plot3d(subs(Sigma=l,soll),t=-10..10,g=-0.05..0,axes=boxed); 




The second solution is: 

>sol2 := subs( 
{rho0=sol2_rho[l],tau=sol2_tau[l],epsilon=sol2_e[l]},expand( 
(f3(t) + f4(t)) 2 ) ): 

>plot3d(subs(Sigma=l,sol2),t=-10..10,g=-0.05..0,axes=boxed); 
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And the unperturbed solution is: 

>sol3 := subs( 
{rho0 2 =-4*g/Sigma,tau=sqrt(-g),epsilon=0},expand( 
(f3(t) + f4(t)) 2 ) ): 

>plot3d(subs(Sigma=l,sol3),t=-20..20,g=-0.05..0,axes=boxed); 
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We can see that the perturbations widen the pulse and reduce its intensity. 
Moreover, the perturbation of the first type splits the pulse. 



We suppose, that the excitation of similar perturbations located 
within the ultrashort pulse dissociates it for the large a, when 
the contribution of the higher-order nonlinear terms is essential 



The important feature of the multiple pulse regimes is the possibility of 
a strong correlation of the pulse parameters in the multipulse complex. This 
is evidence of the pulse interaction. An example of such interaction can be 
illustrated by the following consideration. Let's take the first momentum of the 
Landau- Ginzburg equation by the multiplication to the conjugated field, adding 
the conjugated equation and the consequent integration. The evolution of the 
energy is described by 

>Diff(E(z),z) = 2*g*E(z) + 
int(conjugate(rho(z,t))*Diff(rho(z,t),T(t,2)) + 
rho(z,t)*Diff(conjugate(rho(z,t)),'$'(t, 2)), t=-infinity.. infinity) 
+ 2*Sigma*conjugate(rho(z,t)) 2 *rho(z,t) 2 ; 



E(z) = 2gE(z)- 



( P (Z, t)) (^ P(*. f )) + P( Z > f ) (Qp (P( Z > *))) dt + 2 S M Z ' f )) P( Z > f ) 2 



The second term by the virtue of 

>with(student): 

>intparts(int(rho(z,t)*diff(rho(z,t),'$'(t,2)),t), rho(z,t));# the first term van- 
ishes at infinity 



p(z,t)(^p(z,t))- j{-p{z,t)fdt 



gives: 

>eql4 := Diff(E(z),z) = 2*g*E(z) - 2*int(Diff(conjugate(rho(z,t)),t)* 
Diff(rho(z,t),t), t=-infinity.. infinity) + 2*Sigma*conjugate(rho(z,t)) 2 *rho(z,t) 2 ; 
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eqU := 
^ El : ) = 2 ry E( : ) - 2 / (- WM))) (^ /»(«, *)) * + 2E(^lj) p(*. *) S 



Now let's consider the simplest two-pulse complex: 

>fieldr := 
rhoO* (sech( (t-delta) *tau) +sech( (t+delta) *tau) *cos(phi) ) ;# real part 

>fieldim := rhoO*sech((t+delta)*tau)*sin(phi);# imaginary part, delta is 
the distance, phi is the phase difference 

>print('the spectral term is defined by:'); 

>diff(fieldr,t) 2 +diff(fieldim,t) 2 ; 



fieldr := pO (scch((i -S)t) + sech((t + S) r) cos(<j>)) 



fieldim :— pO sech((i + 6) t) sin(^) 



the spectral term is defined by : 

pO 2 (-scch((t - 5) t) tanh((£ -8)t)t- scch((i + S) r) tanh((t + S)t)t cos(0)) 2 
+ pO 2 sech((i + S) t) 2 tanh((* + 6) r) 2 t 2 sin(0) 2 

The spectral loss for this complex is (second term in eql4): 
>assume(tau>0): 

>s := 2*rho0 2 *tau 2 *( 
int (sech ( (t-delta) *tau) 2 *tanh( (t-delta) *tau) 2 ,t =-infinity. .infinity) 
+ 

int (sech( (t+delta) *tau) 2 *tanh( (t+delta) *tau) 2 ,t=-infinity.. infinity) 
+ 

2*int (sech ( (t-delta) *tau) *tanh ( (t-delta) *tau) *sech( (t+delta) *tau) *tanh 
( (t+delta) *tau) *cos(phi) ,t=-infinity. .infinity)) ; 

>en := int(fieldr 2 + fieldim 2 ,t=-infinity.. infinity); 
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s:=2 P 2 t~ 2 

4 1 4 (ln(%l) e^ 8 T ~ *) - 4 e^ 8 T ~ *) + 6 ln(%l) %1 + 4 + ln(%l)) e< 2 T ~ *) cos(^) 

37 (_ 3e (8r^) +3%1+e (12r-5) _ y T ~ 

%l:=e( 4r ^) 



en := 4 



pO 2 (e^ 4 T ' *) - 1 + e< 2 r ~ *) ln(e< 4 T " *> ) cos(0)) 
T-(e( 4l -'*)-l) 



>plot3d(subs(tau=l/10,s/en),delta=0..40,phi=-Pi..Pi,axes=boxed, 
view=0..0.01,title='spectral loss vs. phase and distance' 




So, there exists the potential well, which can absorb the pulses. 
Also, the last term in eql4 contributes to the interpulse attrac- 
tion due to the loss saturation enhancement produced by the 
pulse merging 



We can conclude, that the analytical treatment allowed the comprehension 
of the basic features of the ultrashort pulse dynamics, though they lie out of 
the quasi-soliton model validity. Such collaboration between numerical and an- 
alytical methods realized by means of the Maple faculties has not only technical 
but also heuristic character. 
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10 Mode locking due to a "slow" saturable ab- 
sorber 

10.1 Analytical theory and linear stability analysis 



In two previous parts we considered the ultrashort pulse formation as result 
of the loss saturation by pulse intensity. This supposes the instant response 
of the saturable absorber on the signal variation. The nonresonant (phase) 
nonlinearity obeys this demand even in femtosecond domain. But the resonant 
nonlinearities are more inertial and the time defining the relaxation of their 
excitation lies in the wide region from 100 femtosecond to milliseconds and 
more. Therefore the interaction of the ultrashort pulse with such structures 
differs essentially from the previously considered. 

Let the pulse duration is shorter then the longitudinal relaxation time. Then 
the loss saturation is caused by pulse energy flux passed through absorber. The 
similar situation had considered for dynamical gain saturation in part 4- Here 
we shall take into consideration simultaneously the dynamical gain and loss 
saturation. These effects can be taken into consideration by the expansion of 



the exponential transmission operator eV "«£/ up to second order on pulse 
energy e (see |2(J]). Here E s is the gain or loss saturation energy. Then the 
basic differential equation is 



>restart: 

with (plots): 

master := diff(rho(z,t),z) = (alpha - g - l)*rho(z,t) - 
chi*alpha*rho(z,t)*int(rho(z,zeta) 2 ,zeta=0..t) + 
alpha*rho(z,t)*(chi*int(rho(z,zeta) 2 ,zeta=0..t)) 2 + 
g*rho(z,t) *int (rho(z,zeta) 2 ,zeta=0. .t) - 
g*rho(z,t)*(int(rho(z,zeta) 2 ,zeta=0..t)) 2 + diff(rho(z,t),t$2) 
delta*diff(rho(z,t),t); 



master := — p(z. t) 
oz 



(a-g-l)p(z,t)-xap(z, t) %1 + a p(z, t) x 2 %1 2 + g p(z, t) %1 



gp( z ,t)%l 2 + (^p(z,t)) + 8(^p(z,t)) 



%i:= / P (z,crdc 

o 
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Here x ls the ratio of the loss saturation energy to gain saturation energy (satu- 
ration parameter), I is the unsaturable loss coefficient, g and a are the saturable 
loss and gain coefficients at pulse peak, respectively, 5 is the pulse delay on the 
cavity round-trip. The form of this equation supposes the normalization of time 
on tf , pulse energy on loss saturation energy. We shall suppose the soliton-like 
form of steady-state solution of master: 

>fl := (t)— >rhoO*sech(t*tau);# soliton form 

f2 := (zeta)— >rhoO*sech(zeta*tau); 
ss := rhs (master): 

subs({rho(z,t)=fl(t),rho(z, zeta) =f2 (zeta) },ss): 
simplify (%): 

expand( numer(%)*2/rho0 ): 
eq := collect( collect( combine(%,trig), 
sinh(2*t*tau)),cosh(2*t*tau) ); 



fl :=t-> pOsech(tr) 

f2 :=C^/50sech(Cr) 

:= (r 2 a-r 2 g + ap0 4 X 2 -t 2 1- g pO 4 + r 4 ) cosh(2 1 r) 

+ (-x a pO 2 t + g pO 2 t - S r 3 ) sinh(2 1 t) - 

ap0 4 X 2 +r 2 a + g pO 4 - t 2 g - t 2 I - 3 r 4 



Since this equation is valid at any moment, we have the system of the algebraic 
equations for the coefficients of hyperbolical functions. 

>eql := coeff(eq,cosh(2*t*tau)); 

eq2 := factor(coeff(eq,sinh(2*t*tau))); 
cq3 := expand(eq-eql*cosh(2*t*tau)-eq2*sinh(2*t*tau)); 



1 2 2 , n4 2 2 i n4 , 4 

eql := t a — t g + a pO \ ~ T <• ~ 9 pO + T 



eq2 := -t (S t z + x a pO 2 - g pO 2 ) 

eq3 := -a pO 4 X 2 + t 2 a + g pO 4 - t 2 g - t 2 I - 3 r 4 
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These equations define the inverse pulse duration r, pulse intensity pO 2 , and 
delay S. Let make some manipulations: 

>eq4 := expand( factor(eql+eq3)/(2*tau 2 ) ); 
eq5 := expand((eql-eq3)/2); 



eq4 '■= —t + a — g — I 



eq5 -.^apO^x -9P0 +2r q 



Hence 



>soll := solve(eq4=0,tau 2 );# solution for tau 
solve(subs(tau 4 =soll 2 ,eq5)=0,rho0 4 ); 
factor (%);# solution for rhoO 4 



soil := a — g — I 



a 2 -2ag-2al + g 2 + 2gl + l 2 
ax 2 ~ 9 



{a-g-lf 
ax 2 - g 



Note that the last solution needs some consideration. The comparison with 
second expression gives for intensity: 



V2(a- g-l) 

\Jg-x 2a 



>sol2 := sqrt(2)*(alpha-g-l)/sqrt(g-chi 2 *alpha):# solution for intensity 

And at last, for delay we have 

>subs( {tau 2 =soll, rho0 2 =sol2},expand(eq2/(2*tau)) ): 
simplify (%): 



143 



numer(%): 

expand(%/ (-alpha+g+1) ) : 
sol3 := solve(%=0, delta) ;# solution for delta 



y/2( X a-g) 

solj := — 

V 3 - " x 2 



The gain coefficient at pulse peak is: 

>assume(tau, positive): 

int(fl(t) 2 ,t=-infinity..O): 
subs({tau=sqrt(soll),rho0 2 =sol2},%):# half-pulse energy 

numer( simplify( alphaO/(l+%) - alpha ) ) = 0;# equation for saturated 
gain, alphaO is nonsaturated gain 



aO v g — a x 2 — ot V 9 — a x 2 — a v2 \J a — g — I = 



>sol := solve (%, alpha): 
tau := 'tau': 

We can plot the pulse intensity and duration versus nonsaturated gain for 
different \- It should be noted, that \ < 1 because for the pulse formation the 
loss saturation has to leave behind the gain saturation. 

>plot3d(subs({l=0.01,g=0.05 

},subs(alpha=sol[l],sol2)),alpha0=0.061..1,chi=0..1, 
axes=boxed,title='pulse intensity vs nonsaturated gain'); 



>plot3d(subs({l=0.01,g=0.05 

},subs(alpha=sol[l],log(l/sqrt(soll)))), 

alpha0=0.061..1,chi=0..1,axes=boxed, 

title= 'logarithm of pulse width vs nonsaturated gain'); 
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logarithm of pulse widthj/s nonsaturated gain 




0.4 
0.6alpha0 



1 1 



>plot3d(subs({l=0.01,g=0.05 

},subs(alpha=sol[l],sol3)),alpha0=0.061..1,chi=0..1,axes=boxed,title= 'log 
arithm of pulse width vs nonsaturated gain'); 



logarithm of pulse widthys nonsaturated gain 




0.4 
0.6alpha0 



1 1 



We can see that the pulse width is decreased by gain growth and decrease of 
X- Additionally there is maximum of the dependence of intensity on saturation 
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parameter \- 

Now we shall analyze the ultrashort pulse stability in framework of linear 
theory (see part 4)- In this case the equation for evolution of exponentially 
growing perturbation £(£) is 

>(alpha-g-l-lambda)*xi(t) - chi*alpha*xi(t)*int(rho(zeta) 2 , 
zeta = 0.. t) - chi*alpha*rho(t)*int(rho(zeta)*xi(zeta), 
zeta = .. t) + 

alpha*xi(t)*chi 2 *int(rho(zeta) 2 ,zeta = .. t) 2 + 
g*xi(t)*int(rho(zeta) 2 ,zeta = .. t) + 
g*rho(t)*int(rho(zeta)*xi(t),zeta = .. t) - 
g*xi(t)*int(rho(zeta) 2 ,zeta = .. t) 2 + 
diff(xi(t),'$'(t,2)) + delta*diff(xi(t),t); 



(a-g-l-X)m- 



X a £(t) %1-xa p{t) / p(Q £(£) d( + a S(t) X 2 %1 2 + 9 £(t) %1 

o 



ui'in I p(omd(- 9 m%i 2 + ^m)+H§- t m) 



%1:= / p(0 2 d( 
o 



Here A is the perturbation's growth increment, its positive value corresponds to 
ultrashort pulse destabilization. An integro-differential character of this equa- 
tion raises the ODE's order therefore we shall use some assumptions relatively 
perturbation's envelope. 

Let consider a long-wave limit for perturbations. In this case the pertur- 
bation's envelope is smooth in compare with pulse envelope. Then we can to 
exclude the integration over perturbation. 

>fl := (t)— >rhoO*sech(t*tau):# soliton form 

f2 := (zeta)— >rhoO*sech(zeta*tau): 
eql := (alpha-g-l-lambda)*xi(t)-chi*alpha*xi(t)*int(rho(zeta) 2 , 
zeta= .. t)-chi*alpha*xiO*rho(t)*int(rho(zeta), 
zeta = ..t)+alpha*xi(t)*chi 2 *int(rho(zeta) 2 , 
zeta = ..t) 2 +g*xi(t)*int(rho(zeta) 2 , 
zeta = ..t)+g*xiO*rho(t)*int(rho(zeta), 
zeta = ..t)-g*xi(t)*int(rho(zeta) 2 , 
zeta = ..t) 2 +difT(xi(t),'$'(t,2))+delta*diff(xi(t),t); 

# xiO is the amplitude of perturbation at t=0 
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eql :=(a-g-l-\)$(t) 



X a £(t) %1 - x a £0 p(t) [ p(C) dC + " £(t) x 2 %1 2 + <?£(*) %1 
Jo 

+ fff0p(*)j[W)dC-5e(*)%l 2 + (^e(t)) + *(^(*)) 

t 

2 



%1:= / p(C) 2 dC 
Jo 

The long-wave approximation allows to neglect the second-order derivation in 
compare with first-order one. 

>value( subs({rho(t)=fl(t), 
rho(zeta)=f2(zeta)},eql-diff(xi(t),'$'(t,2))) ); 



i i vuw xa£{t)p0 2 smh(tT 

(a- g-l- X) £(t) - - 



rcosh(tr) 
X ex £0 pO 2 sech(t r) arctan(sinh(tr)) 
r 
| Q$ft)x 2 p0 4 sinhftr) 2 | g^ft) pO 2 sinh(tr) 
r 2 cosh(iT) 2 tcos1i(£t) 

g £0 pO 2 sech(i r) arctan(sinh(i r)) 
r 
g CWp0 4 sinh(tr) 2 d 

~ r 2 cosher) 2 +<5( ^ C(i)) 

>subs({op(2,%)=-chi*alpha*xi(t)*rho0 2 *tanh(t*tau), 
op(4,%)=chi 2 *alpha*xi(t)*rho0 4 *tanh(t*tau) 2 , 
op(5,%)=g*xi(t)*rho0 2 *tanh(t*tau), 
op(7,%)=-g*xi(t)*rho0 4 *tanh(t*tau) 2 },%); 



(a - 5 - / - A) £(t) - x a £(*) pO 2 tanh(t r) - 

X a £0 pO 2 sech(i r) arctan(sinh(t r)) 
r 
X 2 a £(t) pit tanh(i t) 2 + .9 £(£) pO 2 tanh(t r) + 
g £0 ,o0 2 sech(t r) arctan(sinh(i r)) 



. 9 eWp0 4 tanh(tT) 2 +( 5(^(t)) 
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>dsolve({%=0, xi(0)=xiO} ,xi(t)); 



£0 pO 2 arctan(sinh(u r)) ( «^(«r) -coBh( tt r) )( _ 1/a ^^■, 7 ->") ) 
cosh(ur) 

.coshfur, + smh(uT) , f-, /9 po 2 (%i-a+ X a-apo 2 ) > . . 

— -T- f — r^ — L r f ^ ) (X« - 3 

cosn(w t) 

, -CUT CQ3h(u T) + 8 U T COBh(u t) + I 11 T CQ B h(u x) + A M T CQ 3 h(u T ) + C P Q 4 X 2 3 illh(u t)- 8 P Q 4 3 illh(u t) N 

gV cosh(u t) St > 

1(8 t cosh(w r))du(tanh(t t) - 1) %2 (1 + tanh(i r))^ 1 / 2 ~ (%1 - a + T *°~ gp ° 2) ) 



f 



e t x + g t t + 1 t x + A t x + a pO 4 x 2 tanh(t x)-g pO 4 tanh(f x) 



) + (tanh(tr)-l) %2 £0 



(l + tanh(tr)) ( " 1/2 * 



pO A (%1-g + xa-gpO- 1 ) - 



/ —a t x+g t t+1 t t + A t T+a pO 4 x z tanh(£ x)-g pO a tanh(t x) * 

e ( 3T ) 



/(-1) %2 

%l:=ap0 2 x 2 
1 p0 2 (%l+. 9 - X «-.gp0 2 ) 



%2 := 2 



The last term in this expression has not appropriate asymptotic behavior at 



infinity. As consequence, there are not long-wave excitations in our case 



A short-wave approximation allows to simplify the problem on the basis 
of Riemann-Lebesgue theorem: J^° i(t) e' 7 "'' dt =o(l) ( w~> oo). That is 

the integral J Q p(£) £(£) d( is small value if £(£) is quickly oscillating function 
without steady-state points. In this case (when r<<l ) we have (see |2l|| ): 

>fl := (t)— >rhoO*sech(t*tau):# soliton form 
f2 := (zeta)— >rhoO*sech(zeta*tau): 
eql := 
(alpha-g-l-lambda)*xi(t) - chi*alpha*xi(t)*int(rho(zeta) 2 ,zeta = ..t) + 
alpha*xi(t)*chi 2 *int(rho(zeta) 2 ,zeta = .. t) 2 + 
g*xi(t)*int(rho(zeta) 2 ,zeta = .. t) - 
g*xi(t)*int(rho(zeta) 2 ,zeta= .. t) 2 + 
diff(xi(t),'$'(t,2)) + delta*diff(xi(t),t); 
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eql := (a - g - I - A) £(*) - \ a £(i) %1 + a £(t) x 2 %1 2 + 
.9 £(*) %1 - <?£(*) %1 2 + (^ £(*)) + * (| £(*)) 

%1 := /^C) 2 d( 
Jo 

>subs({rho(t)=fl(t), rho(zeta)=f2(zeta)},eql) ); 

/ , Uf /,N xa^t)pQ 2 smh(tT) 
(a — g — I — A) Hi) -7- — r h 

a £0) X 2 pO 4 sinh(t r) 2 .9 £(*) pO 2 sinh(t r) 
r 2 cosh(tr) 2 rcosh(tr) 

g eWpQ 4 sinhftr) 2 a 2 a 

r 2 coshftr) 2 + fe m + S { m m) 

We can rewrite this expression: 

>(alpha-g-l-lambda)*xi(t)-chi*alpha*xi(t)*rho0 2 * 
tanh(t*tau)/tau+alpha*xi(t)*chi 2 *rho0 4 * 
tanh(t*tau) 2 /tau 2 +g*xi(t)*rho0 2 * 
tanh(t*tau)/tau-g*xi(t)*rho0 4 *tanh(t*tau) 2 /tau 2 + 
diff(xi(t),'$'(t,2))+delta*diff(xi(t),t): 
eq2 := collect (%, tanh); 

Al = expand( coeff(eq2, tanh(t*tau) 2 )/ xi(t)): 
A2 = expand( coeff(eq2, tanh(t*tau))/ xi(t)): 
eq3 := subs( alpha-g-l-lambda=A3,Al*xi(t)*tanh(t*tau) 2 + 
A2*xi(t)*tanh(tau*t) + coeff(eq2,tanh(t*tau),0) ); 



(^^!_Wkg! )tann ( tT) 

T T 

+ (a- g -i-x) m + (^ m) + 5 (| ew) 



e ? 5 := A 1 £(t) tanh(t t) 2 + A2 £(i) tanh(i r) + 45 £(i) + 
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>dsolve(eq3=0, xi(t)); 



_Cihypergeom( 

[%1 + %3 + %2, -%3 + 1 + %2 + %1], [-1 ~ 2r ~ 4%2r + ^ ], 

2 T 

— tanh(tr) + -)(tanh(£r) - 1) %2 (1 + tanh(tr)) %1 + 

hypergcom( 

1 -2%3t + 2t + (5 + 2%1t-2%2t 1 2 %1t + 2%3r + 5 - 2%2t 
[ 2 t ' 2 r J ' 

1 2r-4%2r+j _1 h(ir) + I)( ta nh(ir) - 1) %2 (1 + tanh(tr)) %1 

2 t 2 2 

%1 := RootOf {A3 - A2 + Al +25 _Z r + 4 _Z 2 t 2 ) 

%2 := RootOf (4 _Z 2 t 2 - 2 5 _Z t + A2 + A3 + Al) 

%3 := RootOf (Ai - _Zr 2 + _Z 2 r 2 ) 



>coeffl := allvalues( 
-RootOf(Al-_Z*tau 2 +_Z 2 *tau 2 ) + l+ 

RootOf(A3-A2+Al+2*delta*_Z*tau+4*_Z 2 *tau 2 ) + 
RootOf(4*_Z 2 *tau 2 -2*delta*_Z*tau+A2+A3+Al) ):# 8 coefficients 

coeff2 := allvalues( 
-l/2*(2*RootOf(Al-_Z*tau 2 +_Z 2 *tau 2 )*tau-2*tau+delta- 
2*RootOf(4*_Z 2 *tau 2 -2*delta*_Z*tau+A2+A3+Al)*tau+ 
2*RootOf(A3-A2+Al+2*delta*_Z*tau+4*_Z 2 *tau 2 )*tau)/tau ):# 8 coef- 
ficients 

We have the set of 16 first coefficients of hypergeometric functions, which can 
cause the appropriate asymptotic behavior at infinity if they are equal to neg- 
ative integers (see, in particular, part 4)- Moreover, they are to be the large 
negative integers in order to satisfy to our short-wave approximation. So, we 
investigate the high-level excitations in the "potential well" formed by gain and 
loss saturation. 

>#first coefficients in hypergeometric functions 
c := array(1..16): 

for i from 1 to 8 do 
c[i] := subs({Al = -g*rho0 4 /(tau 2 )+alpha*chi 2 *rho0 4 /(tau 2 ), 
A2 = g*rho0 2 /tau-chi*alpha*rho0 2 /tau, 
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A3 = alpha-g-l-lambda},coeffl[i]): 
od: 
for i from 9 to 16 do 
c[i] := subs({Al = 

-g*rho0 4 /(tau 2 )+alpha*cm 2 *rho0 4 /(tau 2 ), 
A2 = g*rho0 2 /tau-chi*alpha*rho0 2 /tau, 
A3 = alpha-g-l-lambda},coeff2[i-8]): 
od: 

>#16 coefficients produce 16 equations for lambda (N is positive integer) 

s := array(1..16): 

for j from 1 to 16 do 

s[j] := solve(c[j] + N = 0, lambda): 
od: 

># the solutions will be evaluated numerically by variation of chi for differ- 
ent N 

# Attention! This computational block can take a lot of time! 
P := array(1..50,1..16,1..3): 
1 := 0.01: 

g := 0.05: 
alphaO := 0.5: 
Lev := [5, 10, 50]: 
for k from 1 to 3 do 
for j from 1 to 16 do 
for i from 1 to 50 do 
N := Lev[k]: 
chi := i/50: 

alpha := evalf( sol[l] ): 
tau := evalf( sqrt(soll) ): 
rhoO := evalf( sqrt(sol2) ): 

delta := evalf( sol3 ): 
P[i,j,k] := evalf( s[j] ): 
od: 

od: 
print (k); 
od: 

># list of plots 

macro(usercol = COLOR(RGB, 0.8, 5/Lev[k], 0.5)): 
for k from 1 to 3 do 

for m from 1 to 16 do 
p[(k-l)*16+m] := listplot([[n/50,Re(P[n,m,k])] 
$n= 1 . .50] ,color=usercol) : 
od: 
od: 
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> display ({p[ii] $ii=1..40},axes=boxed, 
title='stability increment', view=-2.. 15); 



stability increment 




So, we have only two different solutions for real part of A, when "level's number" 
N is fixed. The negative value of increment corresponds to decaying perturba- 
tions, i. e. stable ultrashort pulse generation. 



Because of the pulse stabilization results from increase of x there 
is problem of the shortest pulse generation (note, that increase 
of x increases pulse width, see corresponding figure above) 



It is of interest that the region of the pulse stability can have an inhomoge- 
neous character (a "kink" on the curves for small N, that is "deeper level" in 
potential well). This fact was confirmed by numerical simulations (see p2|| ) 
and corresponds to destabilization due to excitation of "deeper levels", that is 
the pulse envelope splitting. Perturbations with large N can be interpreted as 
continuous- wave generation of noise spikes. 



10.2 Aberrationless approximation 

Now we shall consider one practical realization of the described here method 
of ultrashort pulse generation. The most attractive sort of "slow" absorber is 
the semiconductor absorber with comparatively short relaxation time and small 
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energy of saturation (see next part). But the absorber, which is based on the 
impure dielectric crystal, is much more usable in the picosecond time domain 
due to its simplicity (as rule, the semiconductor shutter has a very complicate 
structure) and durability. However, the basic disadvantages of the crystalline 
absorbers are the large relaxation time (10 ns - 1 /zs) and large saturation energy 
(up to 1 J / cm 2 ). Here we shall analyze the possibility of the stable ultrashort 
pulse generation in the forsterite solid-state laser with YAG: V 3+ crystalline 
absorber. As the basic method the aberrationless approximation will be used. 
Now we shall express the loss coefficient in absorber as ge"~~^\ that differs 
from phenomenological expression, which was considered in the beginning of 
this part, but is close to it when the expansion on energy is made up to small 
orders of e. Such form of the dependence can be obtained from Bloch's equations 
for two-level absorber (see next part) in noncoherent approximation and in the 
condition of small pulse duration in compare with longitudinal relaxation time 
Ta (Ta = 22 ns for YAG: V 3+ ). The steady-state (there is no dependence on 
z) master equation is (see above): 

>restart: 
master:= 
alpha* (rho(t)-rho(t) *chi*epsilon(t) +rho(t) *chi 2 *epsilon(t) 2 / 
2)-g*(rho(t)-epsilon(t)*rho(t)+rho(t)*epsilon(t) 2 /2) + 
diff (diff (rho(t) ,t) ,t) +delta*diff(rho(t) ,t)-l*rho(t) ; 



master := a (p{t) - p{t) X e(t) + - p(t) X 2 e{t) 2 )- 



g ( P (t) e(t) p(t) + \ p(i) e{tf) + (J^ p{t)) + 5 (^ p{t)) I p{t) 



Let suppose the soliton-like shape of the pulse. Then the evolution of the ultra- 
short pulse parameters can be found as result of the expansion at i-series: 

>fl:=(z,t)— >rhoO(z)*sech(t*tau(z));# pulse amplitude 
f2:=(z,t)- >rhoO(z) 2 *(l+tanh(t*tau(z)))/tau(z);# pulse energy 

lhs_master:=subs({diff(rhoO(z),z)=x,diff(tau(z),z)=y 

},diff(fl(z,t),z)):# left-hand side of dynamical equation 

eq:=collect(series(lhs_master-subs({rho(t)=fl(z,t),epsilon(t)=f2(z,t) 
}, master) , t=0, 3) ,t):#dynamical equation 

eql := coeff(eq,t 2 ): 
eq2 := coeff(eq,t): 

eq3 := coeff(eq,t,0): 
sol := factor(solve({eql, eq2, eq3},{x, y, delta})); 
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fl := (z, t) -> pO(z)scch(ir(z)) 

pO(z) 2 (l + tanh(tT(z))) 



12 := («, t) 



r(*) 



so/ := {3; = — pO(z) 



(2ar(z) 2 - 2apO(zf X T(z) + a pO(z) 4 X 2 -21t(z) 2 - 2 ff r(z) 2 
2gp0(z) 2 r(z) - g pO(z) 4 - 2r(z) 4 ) /r(z) 2 , 
1 ap0(z) 4 x 2 -gp0(z) 4 + 4r(z) 4 



?/ 



2 t(z) 



p0(z) 2 (-g p0(z) 2 -q X t(z) + a p0(z) 2 x 2 + .9 t(z)) , 

" = ^w } 



The equations for the pulse parameters evolution have to be supplemented by 
the equations for the gain ( a^ ct{z)) and the saturable loss (g*-*g(z)) evolution 
at the time intervals >>1/ t ( aO and gmx are the maximal saturable gain and 
loss, respectively, Ta and Tr are the loss and gain relaxation times normalized 
to the cavity period Tcav, Pump is the dimensionless pump, see part 9): 

>eq4:=-4*rhoO(z) 2 *g(z)/tau(z) + (gmx-g(z))/Ta; 

eq5:= 

Pump*(alpha0-alpha(z))-2*chi*alpha(z)*rho0(z) 2 /tau(z)-alpha(z)/Tr; 



, A P0(z) 2 g(2:) gmx-g(z) 

t(z) Ta 

2\a{z) pO(z) 2 a(z) 



eq5 := Pump (aO — ct(z)) — 



t(z) Tr 



Then finally: 



>sys := 

D (g) (z) =eq4,D (a) (z) =eq5,D (rhoO) (z) =subs({alpha=alpha(z) ,g=g(z) } , 

subs(sol,x)),D(tau)(z)=subs({alpha=alpha(z),g=g(z)},subs(sol,y)); 
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# basic systems for evolution of the pulse parameters, gain and loss coeffi- 
cients 



sys := D(g)(z) = -4 



D(o)(z) = Pump (aO — a(z)) — 

D(pO)(z) = ipO(z) 



t(z) Ta 

2 x a(z) p0(z) 2 a(z) 



t(z) Tr 



(2 a(z) t{z) 2 - 2 a(z) pO(z) 2 x t(z) + a{z) P 0{z) 4 x 2 ~ 2 1 t{z) 2 - 
2 g(z) t(z) 2 + 2 g(z) pO(z) 2 t(z) - g(z) pO(z) 4 - 2 r(z) 4 ) /r(z) 2 , 

1 a(z)pO(z)V-g(z)pO(z) 4 + 4T(z) 4 

2 t(z) 

We shall change the saturation parameter %and to search the stationary 
points of the pulse parameter's mapping. These points correspond to the solu- 
tions of 'sys ' with zero left-hand sides 

>st_soll:=solve({eq4,eq5},{g(z),alpha(z)}): 
st_g:=subs(st_soll,g(z)); 

st_a:=subs(st_soll,alpha(z)); 
st_sysl:=[expand(rhs(op(3,[sys]))*2*tau(z) 2 /rho0(z)), 
expand(rhs(op(4,[sys]))*2*tau(z))]: 
st_sys2:={simplify(op(l,st_sysl)+op(2,st_sysl))=0,op(2,st_sysl)=0}; 

st_sys3:=subs({rho0(z)=x,tau(z)=y},subs({g(z)=st_g, 
alpha(z)=st_a},st_sys2)); 



t(z) gmx 
st_g := 



st a :- 



ApO(z) 2 Ta + T{z) 
Pump t(z) Tr aO 



Pump t(z) Tr + 2 \ pO(z) 2 Tr + t(z) 



st_sys2 := 

{2 a(z) t{z) 2 - 2 a(z) p0(z) 2 x t(z) - 

21t(z) 2 - 2g(z)r(z) 2 + 2g(z) p0(z) 2 t(z) - 6t(z) 4 = 0, 

-a(z) P 0{zf x 2 + S(z) p0(z) 4 - 4 r{zf = 0} 
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st_sys3 :— 
Pump y 3 Tr aO 2 Pump y 2 Tr aO x x 



t rj i'\. i o , ,/ : 7'.. i .. }'},.., ■!'.. ■"> ~ , ^,2 



Pump y Tr + 2\x 2 Tr + y Pump y Tr + 2 \ x2 Tr + y 

2 y 3 gmx 2 y 2 jrai x 2 

Ax 2 Ta + y Ax 2 Ta + y 

Pump y Tr aO x 4 % 2 V gmx x A 



Pump y Tr + 2\x 2 Tr + y Ax 2 Ta + y 



4y 4 -0} 



The next procedure will be used for numerical solution of st_ sysS. 

>num_sol := proc(alphaO, gmx, chi,l,Tr,Ta, Pump) 

st_sys := 
{-Pump*y*Tr*alphaO*x 4 *chi 2 / 

(Pump*y*Tr+2*chi*x 2 *Tr+y)+y*gmx*x 4 /(4*x 2 *Ta+y)-4*y 4 = 0, 
2*Pump*y 3 *Tr*alphaO/(Pump*y*Tr+2*chi*x 2 *Tr+y)- 
2*Pump*y 2 *Tr*alphaO*x 2 *chi/(Pump*y*Tr+2*chi*x 2 *Tr+y) 
-2*l*y 2 -2*y 3 *gmx/(4*x 2 *Ta+y)+2*y 2 
*gmx*x 2 /(4*x 2 *Ta+y)-6*y 4 = 0}: 

fsolve(st_sys,{x,y},{x=0..1,y=0..1}): 
end: 

>v := array(1..100): 

for i from 1 to 100 do 
v[i] := num_sol(0.5,0.05,l/(1.36+9*i/100),0.01,300,2.2,0.001) od: 
# the normalization of relaxation times to cavity period is supposed (Tcav 
= 10 ns) 

Now we can plot the logarithm of the pulse duration versus \- 

>with(plots): 

ww:=array(1..100): 
for j from 5 to 100 do ww[j]:=evalf(loglO(l/subs(v[j],y))) od: 
points := {seq([l/(1.36+9*j/100),ww[j]],j=1..100)}: 

fl := 
plot(points,x=0.1..0.6,style=point,symbol=circle,color=red): 

display(fl,TEXT([10,3],"Up=0.0008"),view=2.4..3.6, 
title='Logarithm of pulse duration versus sigma', 
labels= ['saturation parameter',"]); 
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3.6i 



Logarithm of pulse duration versus sigma 




0.2 0.3 0.4 0.5 

saturation parameter 



0.6 



This figure demonstrates the pulse width decrease due to decrease of %, that 
corresponds to the predominance of the loss saturation over gain saturation (see 
previous subsection). For the time normalization to the inverse bandwidth of 
YAG: V 3+ absorption line the pulse duration at x = u -3 is about of 50 ps. 

Now, we shall consider the ultrashort pulse parameters evolution on the 
basis of the obtained system of ODE. This procedure solves the systems by the 
standard operator DEplot 

>with(DEtools): 

ODE_plot := proc(alphaO,gmx,chi,l,Tr,Ta,Pump) 

sys := [D(g)(z) = 
-4*rhoO(z) 2 *g(z)/tau(z) + (gmx-g(z))/Ta, 

D (alpha) (z) =Pump*(alphaO-alpha(z))- 

2*chi*alpha(z)*rhoO(z) 2 /tau(z)-alpha(z)/Tr, 
D(rhoO)(z) = 
l/2*rho0(z)*(-2*l*tau(z) 2 +2*alpha(z)*tau(z) 2 - 
2*alpha(z)*rho0(z) 2 *chi*tau(z)+alpha(z)*rho0(z) 4 *chi 2 - 
2*tau(z) 4 -2*g(z)*tau(z) 2 +2*g(z)*rho0(z) 2 *tau(z)- 
g(z)*rho0(z) 4 )/(tau(z) 2 ), 

D(tau)(z) = 
-l/2*(alpha(z)*rho0(z) 4 *chi 2 -g(z)*rho0(z) 4 + 
4*tau(z) 4 )/tau(z)]: 

DEplot(sys,[rhoO(z),tau(z),g(z),alpha(z)],z=0.. 10000, 
[[rho0(0)=le-5,tau(0) = le-3,g(0)=gmx,alpha(0)=0]], 
stepsize=l,scene=[z,rhoO(z)],axes=FRAME,linecolor=BLACK): 
end: 
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Let vary the saturation parameter \ f° r the fixed pump: 
>display(ODE_plot(0.5,0.05,0.3,0.01,300,2.2,0.001)); 



0.008- 



0.006- 
? 0.004- 



0.002- 

o : - J 
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So, we have the decaying oscillations of the ultrashort pulse amplitude (see part 
7). The decrease of saturation parameter produces 



>display(ODE_plot(0.5,0.05,0.1,0.01,300,2.2,0.001)); 
0.025 : 

0.02- 
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One can see the growth of the ultrashort pulse oscillations (au- 
tooscillations or so-called Q-switch mode locking). This is a 
main destabilizing factor for the considered type of the mode 
locking regime. This result is in agreement with above obtained 
(note, that now we are outside of framework of linear perturba- 
tion theory) 



The important obstacle for the pulse generation by slow absorber is the noise 
growth. In order to investigate pulse stability against noise in framework of the 
considered model we have to add the equation for the evolution of noise energy 
n to sys. 

>ODE_noise := proc(alphaO,gmx,chi,l,Tr,Ta,Pump) 

sys_noise := [D(g)(z) = 
-4*rhoO(z) 2 *g(z)/tau(z) + (gmx-g(z))/Ta, 

D (alpha) (z) = 
Pump*(alpha0-alpha(z))-2*chi*alpha(z)*rho0(z) 2 /tau(z)-alpha(z)/Tr, 

D(rhoO)(z) = 
l/2*rho0(z)*(-2*l*tau(z) 2 +2*alpha(z)*tau(z) 2 - 
2*alpha(z)*rhoO(z) 2 *chi*tau(z)+alpha(z)*rhoO(z) 4 *chi 2 - 
2*tau(z) 4 -2*g(z)*tau(z) 2 +2*g(z)*rho0(z) 2 *tau(z)- 
g(z)*rhoO(z) 4 )/(tau(z) 2 ), 

D(tau)(z) = 
-l/2*(alpha(z)*rho0(z) 4 *chi 2 -g(z)*rho0(z) 4 +4*tau(z) 4 )/tau(z), 

D(n)(z) = 
(alpha(z)-l-(gmx+Ta*(g(z)-gmx)*(l-exp(-l/Ta))))*n(z)]:#see V.L. Kalash- 
nikov et al. Opt. Commun., v.159, 237 (1999) 

DEplot3d(sys_noise,[rhoO(z),tau(z),g(z),alpha(z),n(z)], 
z=0..10000,[[rho0(0) = le-5,tau(0)=le-3,g(0)=gmx,alpha(0)=0,n(0)=le-5]], 
stepsize=l,scene=[z,rhoO(z),n(z)],axes=FRAME,linecolor=BLACK): 

end: >display(ODE_noise(0.5,0.05,0.3,0.01,300,2. 2,0.001), 

title='Noise energy evolution'); 
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Noise energy evolution 
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As one can see from this picture, the noise energy decays, i. e. 
50 ps "auto-stable" pulse is stable against noise, too 



11 Coherent pulses: self-induced transparency in 
the laser 



Here we shall calculate the main characteristics of ultrashort pulses in the solid- 
state laser in the presence of self-focusing in the active crystal and coherent 
interaction with a semiconductor absorber. The good introduction in this prob- 
lem and the numerical simulations of the nonlinear dynamics in the laser with 
coherent absorber can be found in the article |23| . 

We shall suppose, that as saturable absorber the semiconductor structure is 
used (like GaAs/AlGaAs quantum- well structures). The key characteristics of 
this absorber, which correspond to the experimental dates, are the energy flux of 
the loss saturation E a = 5* lO^ -5 ) J/ cm 2 , and the time of coherency t co h = 50 

fs. These parameters define the parameter q = J 2 ^tcoh ( n ^ s tne refractivity 
index, c is the light velocity in the vacuum, eO is the dielectric constant). 

>restart: 

with (plots): 

with(DEtools): 

parameter_q := proc() 
local epsilonO,n,c,Ea,q,tcoh: 
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epsilonO := 8.85e-14: 
n := 3.32: 
c := 3el0: 
Ea := 5e-5: 
tcoh := 5e-14: 

sqrt(epsilonO*n*c/(2*Ea*tcoh)): 
end: 

q=evalf(parameter_q());# The dimension of q is [cm/V/s] 

q = .4198714089 10 8 

The corresponding dipole moment is d = q*h/(2* tt): 

>d=evalf(subs({h=6.63e-34,q=parameter_q()},q*h/(2*Pi)));# The dimen- 
sion of d is [coulomb* cm] 

d/e=evalf(rhs(%)/1.6e-19*le7);#heree is the elementary charge, the dimen- 
sion of d/e is [nm] 

d= .4430471653 1CT 26 

d 

- = .2769044783 
e 

For a given E a and the photon energy 2.5* 10^ 19 ^ J (the wavelength is equal 
to 800 nm) the absorption cross-section S = 5* 10^~ 15 ^ cm 2 . The free carrier 
density TV" = 10 19 cm^~ 3 ^ in the semiconductor produces the loss coefficient 
gam_ abs = 0.05 for the length of the absorber's layer z_ abs = 10 nm. Other 
important characteristic is p = 2 tt N d q uj/c = ^ 7r 2 N d 2 lo/(c h) ( u> is the 
field carrier frequency). These parameters are connected with a loss coefficient: 
gam_ abs/z_ abs = NY, = p *tcoh, so 

>p=evalf(subs({tcoh=5e-14,N=lel9,Sigma=5e-15},N*Sigma/tcoh)); 
#The dimension of p is [1/cm/s] 

p= .1000000000 10 19 

As the laser we shall consider Ti: sapphire solid-state laser that is the typical 
laser for femtosecond generation. Its inverse gain bandwidth defining the min- 
imal pulse duration tf = 2.5 fs. We shall use the next normalizations: 1) the 
time is normalized to tf, 2) the field is normalized to 1/q/ tf = 0.95* 10 7 V/cm; 
3) the field intensity is normalized to eon c/(2q tf) 2 , that is 

>intensity_normalization_parameter := procQ 
local epsilonO, n,c,q,tg, par: 
epsilonO := 8.85e-14: 
n := 3.32: 
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c := 3el0: 
tf := 2.5e-15: 
q := parameter_q(): 

epsilon0*n*c/(2*q*tf) 2 : 
end: 
evalf(intensity_normalization_parameter());# The dimension is [W/cm 2 ] 

.2000000000 10 12 

For such normalization the parameter of the gain saturation for Ti: sapphire 
(the energy of the gain saturation is 0.8 J/cm"2) is £ = 6.25* 10*- -4 '. As it 
will be seen, the principal factor in our model is the inverse intensity of the 
saturation of diffraction loss in the laser (so-called Kerr - lens mode locking 
parameter) a = 0.14 (parameter of the diffraction loss saturation is 10 7 W and 
the mode cross section of the Gaussian beam is 30 //m). The corresponding 
SPM coefficient in 1 mm active crystal is j3 = 0.26. 

The response of the semiconductor absorber is described on the basis of 
Bloch equations [£4| . In the absence of the field phase modulation and tuning 
from the medium resonance frequency and in the absence of the relaxation over 
time interval, which is comparable with the pulse duration, the evolution of the 
polarization (components a(t) and b(t)) and the population difference between 
excited and ground states w(t) obey to the following equations: 

>bloch_l := diff(b(t),t)=q*rho(t)*w(t); 

bloch_2 := diff(a(t),t)=0; 
bloch_3 := diff(w(t),t)=-q*rho(t)*b(t); 



d 
block _1 := -q-b(i) =qp(t)w(t) 




d 
bloch 2 := — ait) = 
dt 




d 

bloch 3 := — wit) = -qpit) bit) 
~ ot 



The solutions of this system are 

>sol_bloch := 

dsolve({bloch_l,bloch_2,bloch_3,w(0)=-l,b(0)=0,a(0)=0}, 
{a(t),b(t),w(t)}): 
sol_a := subs(sol_bloch,a(t)); 

sol_b := subs(sol_bloch,b(t)); 
sol_w := subs(sol_bloch,w(t)); 
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sol _a := 

sol _b := — sin(g / p{u) du) 
Jo 

sol_w := — cos(q / p{u) du) 
Jo 



The argument of sin/cos is the pulse area ip: 

>b(t):= -sin(psi(t)); 
a(t):= 0; 



b(i) := -sm(ip(t)) 

a(i) := 



The linear distributed response of the laser with the gain coefficient a, loss 
coefficient gam, frequency filter with the inverse bandwidth tf, laser dispersion 
element with dispersion coefficient hi is described by the terms in the right hand 
side of the wave equation (see Parts 3, 5): 

>Laser_linear := 

alpha*rho(z,t)-gam*rho(z,t)+tf 2 *diff(rho(z,t),t,t) + 
I*k_2*diff(rho(z,t),t,t); 



Laser linear 



d 2 d 2 

ap(z, t) - gamp(z, t) + tf 2 (— p{z, t)) + Ik_2{—p{z, t)) 



The response of nonlinear laser factors is 

>Laser_nonlinear := 

sigma*rho(z,t)*conjugate(rho(z,t))*rho(z,t)- 

I*beta*rho(z,t)*conjugate(rho(z,t))*rho(z,t); 



Laser _nonlinear := a p(z, t) (p(z, t)) — I (3 p(z, t) (p(z, t)) 



Here a and gam are the dimensionless values, that supposes the normalization 
of length z to the length of the optical medium (active medium in our case) , gam 
includes not only scattering loss into optical elements, but also the output loss 
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on the laser mirror, a and (3 have dimension of the inverse intensity, i. e. \p\ 
is the field intensity (we does not write the factor e0*n*c/2, which corresponds 
to the transition field 2 -> intensity). 

As the result, the master integro-differential equation for the field evolution 

is 

>master_l := diff(rho(z,t),z)+diff(rho(z,t),t)/c = 
subs(N=N*z_abs,-2*Pi*N*d*diff(a(t),t)/c+2*Pi*N*d%(t)*omega/c) + 
Laser_linear+Laser_nonlinear;#see Part 3 



,d , .. 4fP{z,t) „irNz abs d sin libit)) to 
master _1 := (— p(z, t)) + ^^ - -° = — 



d 



z c 



9 ty 

ap(z,t)- gamp{z,t) + tf (—^p(z,t)) + 



Ik -% (^2 P& *)) + Q P^ f ) 2 M*' *)) ~ J M^ f ) 2 (/"(*> *)) 



Let transit to the differential equation: 

>assume(q,real): 

master_2 := expand(subs(rho(z,t)=diff(psi(t),t)/q,master_l)); 



■Srjiplt) ir N z absdsm(ib(t))w 
master 2 := dt J = -2 = vyv ;; 



q c c 

777' '-;) gam(^^(t)) | tf (jg 



«(!#)) gam am) , tf 2 {&*m) 



i k _2(§^m) . g(-§- t m?(i t m) , -iHi t m?{§- t m) 



We shall consider only steady-state field states, i. e. -^ p(z, t) = 0, and to 
introduce the time delay on the cavity round-trip 5 

>master_3 := rhs(master_2)+delta*diff(psi(t)/q,T(t,2)); 
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■n N z abs dsm.(ip(t)) u> 

master 3 := —2 h 

c 

aj^m) gam(%m) ] */ 2 (|^(*)) 

(f (f q~ 



q q d 



t0(mW)r(mW)) , <*(&iH*)) 



<f 3 



Let introduce the parameter A, which can be 1) ratio of mode cross-section on 
active medium to one on semiconductor absorber, or 2) coefficient of refractivity 
(for field amplitude) of multilayer mirror on the semiconductor absorber (so- 
called semiconductor saturable absorber mirror - SESAM). Then 

>master_4 := 
subs({sigma=sigma/lambda 2 ,beta=beta/lambda 2 },master_3); 



7r N z abs dsw.(ilj(t))u> 

master _4 := -2 = — — h 

c 



Ik_2{&m) , <r{&W)) 2 (mW)) 



q~ A 2 q~ 3 

+ 



-i/3(& m?{§- t m) , *(^^(*)) 



2 „~3 



A 2 9 



Now let transit to coordinates 'pulse amplitude - pulse area' and eliminate fast 
saturable absorber, GDD and SPM. These suppositions lead to the second-order 
ODE. That is 

>master_5 := 

collect (expand(subs({k_2=0,beta=0,sigma=0 

},master_4)/(-2*Pi*N*d*omega*tf*z_abs/c)), 

cliff (psi(t) ,t ) ) ; # reduced master equation coefficients 

subl := al=-l/2*c/(Pi*N*d*omega*tf*q*z_abs): 

sub2 := a2=-l/2*c*delta/(Pi*N*d*omega*tf*q*z_abs): 
sub3 := a3=-l/2*c*alpha/(Pi*N*d*omega*tf*q*z_abs) + l/2*c*gam/ 
(Pi*N*d*omega*tf*q*z_abs): 
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master_6 := al*diff(psi(t),'$'(t,3))+a2*diff(psi(t),'$'(t,2)) + 
a3*diff(psi(t),t)+sin(psi(t)); 

dsolve(master_6=0,psi(t));# try to solve master_6 and find a very useful 
change of the variables 

master_7 := 

expand(numer (simplify (subs({al=rhs(subl),a2=rhs(sub2),a3=rhs(sub3) }, 

al*(diff(rho(psi),psi,psi)*rho(psi) 2 +diff(rho(psi),psi) 2 * 

rho(psi))+a2*diff(rho(psi),psi)*rho(psi)+a3*rho(psi)+sin(psi))))/(-c)); 

# use the founded change to reduce the order of ODE 



master 5 :- 



ca 



— + 



-cgam 



2 it N duj tf z_abs q~ tt N du tf z_abs q~ I dt 



im) 



sin(^(t)) 1 t/c(^#)) 1 cS(^m) 



tf 



2 tt N duj z_abs q~ 2 tt N duj tf z _abs q~ 



master _6 := al {& V(*)) + a % (& ^(*)) + a3 (& ^(t)) + sin(^(t)) 



?/;(£) = _a&where 



v 9 a 



,_b(_ a ))_b(_ a ) 2 + — 



((- _b(_a)) 2 _b(_a) aJ + a2 (- _b(_o)) _b(_a) + a5 _b(_a) 



sin(_a)) = 



{_M_a) = ^W, _o = ^(*)}, |t = y h } a) d_a + _Cl,i/}(t) = _a 



master _ 7 := (— p(V>)) p(V>) 2 + (^ p(V0) 2 /°W>) + <$ (^ p(V0) p(V0+ 



p(V0 a: — p(V') ^goti — 



2 sin(^) tt N dio tf q~ z_abs 



As it is known (see, for example, [Q) the propagation of the extremely short 

laser pulse in the coherent absorber causes the effect of the self-induced trans- 
parency, when the pulse does not suffer the decay and there is not a transforma- 
tion of the pulse shape in the absorber. This effect is the result of the coherent 
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interaction of the pulse with the atoms and is described on the basis of Bloch 
equations. In the beginning we shall consider the steady-state ultra-short pulse 
in the presence of the coherent interaction with absorber, but without any las- 
ing factors. In this case, the modified master equation master _ 7 contains only 
two terms: the term describing the pulse time delay and the term correspond- 
ing to the coherent polarization response. Note, that p is the field amplitude 
multiplied by parameter q, ip is the pulse area. 

>odel := delta*rho(psi)*diff(rho(psi),psi) - p*sin(psi); 



odel := 5 (— p(ip)) p(ip) -psm(tf>) 



The natural initial condition is p(0) = 0, that is before pulse front its amplitude 
and area are 0. Then solutions of the master equation are 

>sol := dsolve({odel,rho(0)=0},rho(psi)); 



sol := 



pW = 



^/S(-2pcos(tP) + 2p) 



pW = - 



y/5(-2pcos(ip) + 2p) 



As one can see, only one solution ( S > 0) is physical, because of the amplitude 
is to be real. 

The distinct pulse velocity defines the pulse amplitude. This fact is illus- 
trated by the next figure: 

>animate(evalf(abs(subs(p=5e-4,subs(sol[l],rho(psi)))/10)), 

psi=0..2*Pi,delta=0.3e-5..1e-4, 

frames=100,axes=boxed,color=red, 

labels= ['pulse area','y, MV/cm'],title='Pulse amplitude versus its area'); 
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Pulse amplitude versus its area 




2 3 4 

pulse area 



The figure demonstrates, that the full pulse area is equal to 2 
it. Such pulse is named as 2 n - pulse or 2 n - soliton and 
the process of its formation is the so-called phenomenon of the 
self-induced transparency 



More natural representation of the coherent solitons in our case is produced by 
the transition to the coordinates "pulse area - time". Then the master equation 
(see master _4) can be written as follows: 

>ode2 := diff(psi(t),t$2)=a*sin(psi(t));# analog of nonlinear Klein-Gordon 
equation 



ode2 := —, r ijj(t) = asin(ip(t)) 
ot z 



Here a=p/ S. This is the equation, which is analog of the equation of the 
pendulum rotation (the angle variable x in our case is the ultra-short pulse 
area) and the angle is measured beginning from the upper equilibrium point if 
S > (see above). 

>ode2 := diff(x(t),t$2)=a*sin(x(t)); 

ode2 := Jp-x(i) = asin(x(t)) 
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The solution of this equation is well known. 
>soll := dsolve(ode2,x(t)); 



soil := / . d a-t- C2 = 0, 

J v /-2acos(_a) + CI ~ 

d a-t- C2 = 



A/-2acos(_a) + _C1 



Make an explicit integration: 

>sol2_a := value(lhs(soll[l])); 
sol2_b := value(lhs(soll[2])); 



sol2_a := 2^/-(-±a%l + 2 a + _C1 ) (-1 + %1) \J\ - %1 



-4a%l + 2a+ O/ . . „, ,1 



EllipticF(cos(-x(i)), 2./ — ) 



2a+_Cl ' v 2 V 2a+ CJ 

(4acos(ix(i)) 4 -6a%l + 2a- _C1%1 + _C1 
,1 



sin(- x(t)) V-4a%l + 2a+_C*i) - t - _C2 
%l:-cos(ix(i)) 2 



'l-cos(ix(/i)^ 
W2 6 := -2 ' 



\ 



-4 a cos(- x(t)) 2 + 2 a + _ CI 

\ , n %1 ) 

2 a + Gi 



sin(- x(i)) W-4 a cos(- x(t)) 2 + 2 a + _ CI 



-t- _C2 
%1 := EllipticF(cos(-x(i)), 2 



2 w " y2a+_Ci 

The result is expressed through elliptic integrals. The simplification of the 
radicals produces: 

>sol3_a := simplify(sol2_a+t+_C2, radical, symbolic)=t+_C2; 
sol3_b := simplify(sol2_b+t+_C2, radical, symbolic)=t+_C2; 
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solS _a := 

^ l-co S (lx(t))^EllipticF(cos(lx(t)), 2 ^± -^ -) 
- 1 



v /2a + _C*isin(-x(t)) 

= t+ C2 



solS b 



1 - cos(l X (tW EllipticF(cos(l x(f )), 2 ^ ° ± ~ -g^ 
,2 2 zo+ Ci 



.1 



V / 2a+_Cisin(-x(t)) 



So, we have: 

>sol3_a := 

2*EllipticF(cos(l/2*x(t)),2*sqrt(a*(2*a-_Cl))/(2*a-_Cl))/ 
sqrt(-2*a+_Cl) = t+_C2; 

sol3_b := 
-2*EllipticF(cos(l/2*x(t)),2*sqrt(a*(2*a-_Cl))/(2*a-_Cl))/ 
sqrt(-2*a+_Cl) = t+_C2; 



EllipticF(cos(lxft)),2 ^ a 2 (2 ;-- gJ) ) 

solS a := 2 — — = = t + C2 

v /-2a+_Ci 

EllipticF(co S (lx(t)),2 ^ a 2 (2 ^- gi) ) 

solS b := -2 =- — = =t+ C2 

v / -2a+ _C1 

Now we define the initial conditions. Let suppose, that x(0) = n, dx(0)/dt 
= a, where a is the some positive value (we measure the time from lower 
equilibrium point). Then 

>i_Cl := 

solve(simplify(subs({dirT(x(t),t)=alpha,x(t)=Pi}, 
expand(simplify(diff(lhs(sol3_a),t))))) = l,_Cl); 

i_Cl := 
solve(simplify(subs({diff(x(t),t)=alpha,x(t)=Pi}, 
expand(simplify(diff(lhs(sol3_b),t))))) = l,_Cl); 
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CI :=2a + a 2 
CI :=2a + a 2 



The second constant of integration can be found as follows: 

>i_C2 := simplify(subs({_Cl=i_Cl,x(t)=Pi},lhs(sol3_a)))=C2; 
i_C2 := simplify(subs({_Cl=i_Cl,x(t)=Pi},lhs(sol3_b)))=C2; 

i_C2 :=0= C2 
i C2:=0=C2 



The implicit solution is 

>sol4_a := 

simplify (subs(_Cl=i_Cl,lhs(sol3_a)), radical, symbolic) =t+lhs(i_C2); 
sol4_b := simplify(subs(_Cl=i_Cl,lhs(sol3_b)), radical, symbolic) 
=t+lhs(i_C2); 

EllipticF(cos(- x(t)), 2 — -) 

sol4_a :=2 * ^ = i 

a 

EllipticF(cos(-x(i)), 2^—) 
sol4_b:= -2 * S = i 



Let consider a special situation, when 2 ¥-§-= 1: 

>sol5_a := 

solve(simplify(lhs(subs(2*sqrt(-a)/alpha=l,sol4_a)))=rhs(sol4_a),x(t)); 
sol5_b := solve(simplify(lhs(subs(2*sqrt(-a)/alpha=l,sol4_b))) 
=rhs(sol4_b),x(t)); 

sol5 _a := 2 arccos(tanh(— ta)) 
sol5 _b :— 2ir — 2arccos(tanh(- ta)) 



>if sign(sol5_a) < then # we choose only the root corresponding to grow- 
ing area 

sol := sol5_a 

else 
sol := sol5_b 
fi: 
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The dependence of the pulse area on the time is: 

>pulse_area := 
subs(alpha=solve(2*sqrt(p/delta)/alpha=l, alpha), sol); 

animate(evalf(subs(p=5e-4,pulse_area)), 
t=-10..10,delta=0.3e-5..1e-l, 

frames=100,axes=boxed,color=red,labels=['time, t/tf ,'psi'], 
title='Pulse area versus time'); 
#time is normalized to tf 



P 
pulse_area := 2n — 2 arccos(tanh(£ , —)) 

V 



Pulse area versus time 




-2 2 
time, t/tf 



8 10 



The pulse profile is: 

>field := value(simplify(convert(diff(pulse_area,t),sincos))); 



/ 



csgn 



cosh 



t 2 p 



\\ 



field := 2 ■ 



YVU 



cosh(t w-) 
V o 
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This is the coherent soliton with duration 



Zp — 



and amplitude 



2 



>animate(evalf(subs(p=5e-4,field)/10),t=-10..10, 
delta=0.3e-3..1e-l,frames=100,axes=boxed,color=red, 
labels=['time, t/tf','rho,MV/cm'],title='Pulse envelope'); 

Pulse envelope 




We have shown that there is the coherent soliton with sech - 
shape profile in the condition of the coherent propagation in the 
media described on the basis of Bloch equations 



This pulse may be described in the coordinates 'field - area' or 'field - time'. 
The first is formally very simple, but the second representation is physically 
more obvious and corresponds to the model of nonlinear pendulum. 

Now we return to the laser model (master _ 7). This second-order nonlinear 
nonautonomous ODE can be solved numerically (pz = 2n Nd 2 (uj *tf) *z_ abs/(c 
h) = gam_abs* j^, we supposed gam_abs = 0.01) 

>de := subs({ 
alpha=0.1, 
gam=0.04, 

delta=0.0042, 
pz=5e-4}, 
subs({op(6,master_7)=-pz*sin(psi),rho(psi)=rho(psi)/10},master_7)): 
#here 10 results from the time normalization to tf (i.e. the field is measured 
in l/(q*tf) [MV/cm] - units) 
fig:= 
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DEplot([de=0],rho(psi),psi=0.01..1.985*Pi,[[rho(Pi)=0.76*10, 
D(rho)(Pi)=le-15]],rho=0..0.76*10,stepsize=0.001,linecolor=green) 
display (fig, labels= ['pulse area','rho, MV/cm'], 
title='Pulse amplitude versus its area',view=0..7.6); 



Pulse amplitude versus its area 



E 
> 




3 4 

pulse area 



Let compare this result with the secft-shape profile (blue color) corresponding 
to the coherent soliton: 

>plot(subs(aml=0J6,aml*sin(psi/2)*10),psi=0..2*Pi,color=blue): 
display(fig,%,labels=['pulse area','rho, MV/cm'], 
title='Pulse amplitude versus its area'); 
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Pulse amplitude versus its area 



rho, MV/cm4- 




We obtained 2 it - pulse, but, as it can see from previous figure, 
such pulse is not quasi-soliton with sec/i-shape profile (which is 
shown by the blue color, see below about connection between 
coordinates 'field - area' and 'field - time' in this case) 



The field 10 6 V/cm corresponds to the intensity in vacuum 1.3 GW/ cm 2 , that 
is the typical intracavity intensity for the mode-locked solid-state laser. Now 
make a try for obtaining of an approximate solution of the master equation. We 
shall use the solution form, which is typical for the analysis of the equations 
describing the autooscillations (harmonic approximation) : 

>approx_sol := 
aml*sin(psi/2)+am2*sin(psi);# approximation 

f7 := numer(combine(expand(subs(rho(psi) = 
approx_sol,master_7)),trig)); 
# substitution into initial equation 

approx _sol :— ami sin(- tp) + am2 sin(i/>) 
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/7 := 4 5 ami 2 sin(V>) c + 8 8 am2 2 sin(2 ?/>) c + 12 5 amJ am^ sin(- tp) c 
— 4 S ami am2 sin(- ?/>) c + 11 ami am2 sin(2 i/>) c — 16 ga?7i am2 sin(^) c 

— 16 gam ami sin(— ip) c + 16 a am2 sin(ip) c + 16 a ami sin(— ip) c 

3 5 

— 9 ami am2 sin(— ip) c + 17 ami am2 sm(—ip)c — 

32 sin(-0) ir N duj tf q~ z _abs — 8 am2 sin(^) c — 10 ami am2 sin(-0) c + 

3 11 

2 ami sin(— ip) c — 2 ami sin(— ^) c— 10 ami sin(— ip) am2 c + 

8 am2 i sin(3 tp) c 



We have to collect the coefficients of sin( ip/2) and sin( tp)'- 



>f8 := coeff(f7,sin(psi/2)); 
f9 := coeff(f7,sin(psi)); 



/* 



ami = 2 \/W 

V o 



—4 5 ami am2 c — 16 gam, ami c + 16 a ami c — 2 ami c — 10 ami am2 c 

f9 := 4 5 a?ni c — 16 garn ami? c + 16 a ami? c — 32 7r N dui tf q~ z _abs 
-8 am2 i c - 10 ami 2 am2 c 

Note, that in the absence of the lasing factors approx_ sol is the exact solution 
with the parameters am2 = 0, 



Below we will suppose, that 

the approximate solution is close to the symmetrical shape, i. e. am2 = 0. 
Then 

>fl0 := expand(factor(subs(am2=0,f8))/(-2*c*aml)); 
fll := subs(am2=0,f9); 
symmetrical_sol_l := 
allvalues(solve({flO=0,fll=0},{aml,delta})); 

flO := ami + 8 gam — 8a 
fll := 4<5 ami 2 c - 32 tt N dto tf q~ z _abs 



symmetrical _sol _1 := { 



s = 


7r N du tf q~ z abs 


ami - 




= \J— 8 gam + 8 a 


c (—gam, + a) 



, irNdu)tfq~z abs 

{6 = -^ — , ami 

c (—gam + a) 



■\f— 8 gam + 8 a} 



176 



So, the pulse amplitude is defined by the laser parameters, al- 
though the relation between the pulse amplitude, duration and 
delay corresponds to coherent soliton (see above) 



Now return to the coordinates 'amplitude - time' for the harmonic approxima- 
tion. 

>symmetrical_sol_2 := 

dsolve(diff(psi(t),t)-subs({psi=psi(t),am2=0}, 
approx_sol) ,psi(t)) ; 



symmetrical _sol _2 := 

/It aml+l/2 CI ami) _, 

1 -I- P^ ami + _ CI ami ) ' 1 _i_ p(t aml-\-_Cl ami) 



(l/2t aml+l/2 _C1 ami) _ (t aml+_Cl ami) , i 
ljj{t) = 2arctan(2 ,/ taml + CI ami) ' I , t >(taml+ CI ami) ' 



The normalized field envelope is: 

>symmetrical_sol_3:= 

simplify (diff(subs(symmetrical_sol_2,psi(t)),t)); 

ami e (i/2 ami (t+_Ci)) 
symmetrical sol 3 := 2 -. —7—, — 777TT 

y — — e {aml (t+_Cl)) _j_ I 



The initial condition is p(0) = ami. Then 

>in_C := solve(subs(t=0,symmetrical_sol_3)=aml,_Cl); 

m_C :=0, 

And finally we have: 

>symmetrical_sol := subs(_Cl=0,symmetrical_sol_3); 



„ ami e W2tami) 
symmetrical sol := 2 7- -r 

a — g(t ami) 1 ^ 



Naturally, this is a sech - shape pulse with the amplitude 
the duration 



ami = ^ a . 7 

qtj 



and 



tp 



19 V2("-7) 



The pulse profile in the dependence on 



the gain coefficient can be shown by the next function: 
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>plot3d(subs(gam=0.04,2*sqrt (2* (alpha-gam))* 

sech(t*sqrt(2*(alpha-gam))/2.5e-15)*10),t=-le-13..1e-13, 

alpha=0.04..0.1,axes=boxed,orientation=[290,70],labels= 

['t, s', 'alpha', 'rho'],title='Pulse amplitude (MV/cm) versus time'); 

Now we can found the dependence of the pulse parameters on the critical laser 
parameter, that is the pump. For this aim, we have to express the gain coefficient 
from the pump intensity. We shall suppose, that active medium operates as a 
four level scheme. In this case the steady-state saturated gain is described as 
follows ||: 

>alpha=Pump*alphamx/(Pump+tau*Energy+l/Tr); 

Pump alphamx 



a 



Pump + t Energy -\ 

Tr 



Here Pump = a _ab*Tc*Ip/h* v is the normalized pump intensity, a _ab is 
the absorption cross-section at the pump wavelength, Tc is the cavity period, 
Ip is the pump intensity, h* v is the pump photon energy, alphamx is the 
maximal gain, Energy is the normalized pulse energy, Tr is the gain recovery 
time normalized to Tc (dimensionless Tr = 300 for Ti: sapphire laser with 
cavity period 10 ns). 

For harmonical approximation: 

>Energy=2*aml 2 *tp:#pulse energy 

fl2 := Pump*alphamx/(Pump+tau*Energy+l/Tr)-alpha: 
fl3 := 
numer(simplify(subs( 

{aml=2*sqrt(2*(alpha-gam)),tp=l/sqrt(2*(alpha-gam)) 
},subs(Energy=2*aml 2 *tp,fl2)))): 
alpha_sol := solve(fl3=0, alpha): ^solution for the saturated gain 

The dependence of the pulse duration (two physical solutions correspond to two 
different pulse energy) versus dimensionless pump coefficient is 

>fig := plot({ 

Re (evalf (subs (lambda= 1 ,subs ( 

{alphamx=0.1,Tr=300,tau=6.25e-4/lambda 2 ,gam=0.01 
},subs(alpha=alpha_sol[l], 2. 5/sqrt (2* (alpha-gam))))))), 
Re (evalf (subs (lambda= 1 ,subs ( 

{alphamx=0.1,Tr=300,tau=6.25e-4/lambda 2 ,gam=0.01 
},subs(alpha=alpha_sol[3], 2. 5/sqrt (2* (alpha-gam)))))))}, 
Pump=0. 0005. .0. 005, axes=boxed,labels=['Pump, a.u.', 'tp, 
fs'],title='Pulse duration versus pump',color=red): 
display (fig, view=5.. 40); 
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Pulse duration versus pump 




0.001 



0.002 0.003 
Pump, a.u. 



0.004 



0.005 



As we can see, the pump growth decreases the pulse duration 
up to sub-10 fs region 



The quasi-soliton in the absence of the coherent absorber in the laser with self- 
focusing is described in parts VII, VIII. Now we shall demonstrate, that the 
essential features of the lasing in the presence of both factors is the possibility 
of the quasi-soliton generation (compare with above discussed situation). We 
shall search such solutions. 

>assume(t,real): 
assume(tp,real): 

aO := 2/tp: # This is the pulse amplitude 
sol := int(a0/cosh(t/tp),t): # This is the pulse area psi 
subs({disp=0,beta=0,psi(t)=sol,k_2=0,beta=0},master_4): 
expand (%): 
numer(%): 
eql := expand(%/(4*exp(t/tp))); 
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eql := -2irN z_absdutp~ 3 q~ 3 A 2 + 

2tt N z_absduj tp~ 3 q~ 3 X 2 (e (l H) 4 + actp~ 2 q~ 2 A 2 

+ 2actp~ 2 q~ 2 X 2 (e { ^" ) ) 2 + qc £p~ 2 g~ 2 A 2 (e^) 4 - gam ctp~ 2 q~ 2 A 2 

-2gamctp~ 2 q~ 2 A 2 (e (l H) 2 - gamctp~ 2 q~ 2 X 2 (e^) 4 + tf 2 cq~ 2 X 2 

~6tf 2 cq~ 2 X 2 (e { W ) f + tf 2 (e^)4 cq ~ 2 X 2 + 16<j(e { W ) )2 c + 

5 c tp~ q~ 2 X 2 -5ctp~ q~ 2 X 2 (e^) 4 



Collect the terms with equal degrees of exp(t/tp). As result we obtain the 
equations for the pulse and system parameters. 



>el := coeff(eql,exp(t/tp) 4 ); 

e2 := coeff(eql,exp(t/tp) 2 ); 
e3 := expand(eql-el*exp(t/tp) 4 -e2*exp(t/tp) 2 ); 



el := 
2ir N z _abs duj tp~ 3 q~ 3 X 2 +actp~ 2 q~ 2 X 2 - gamctp~ 2 q~ 2 X 2 + tf 2 cq~ 2 X 2 

-5ctp~ q~ 2 X 2 



e2 := 2actp~ 2 q~ 2 X 2 - 2 gamctp~ 2 q~ 2 X 2 -6tf 2 cq~ 2 X 2 + 16 ere 

e3 := 
-2itN z_absduj tp~ 3 q~ 3 X 2 + a c tp~ 2 q~ 2 X 2 - gam c tp~ 2 q~ 2 X 2 + tf 2 c q~ 2 X 2 

+ Sctp~ q~ 2 X 2 



>e4 := simplify (el-e3); 
e5 := simplify (el+e3); 
e6 := simplify (e2-e5); 

e4 := 4nN z _abs duj tp~ 3 q~ 3 X 2 -25ctp~ q~ 2 X 2 

e5 := 2actp~ 2 q~ 2 X 2 -2gamctp~ 2 q~ 2 X 2 + 2 tf 2 c q~ 2 X 2 

e6 := -8 tf 2 c q~ 2 X 2 + 16 a c 



>allvalues(solve({e4=0, e5=0,e6=0},{tp, delta, sigma})); 
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tp' 



1 



—gam + a 



tf,s 



tf 7r N z_abs dui , 
c (—gam + a) 



1 



2 „~2 \2 



tfq~ Z \ 



{tp~ 



1 



-^aTTi 



tf, 6-. 



tf it N z _abs du> q~ 
c (—gam + a) 



- = ^/V 2 A 2 } 



We see the essential differences from the previous situation: 1) 
there is the pulse with sec/i-shape (quasi-soliton); 2) the pulse 
exists, when a < 7, i. e. the linear loss exceeds the saturated 
gain. This is an essential demand to the pulse stabilization and 
breaks the limitations for the loss coefficient in the semiconduc- 
tor absorber; 3) the quasi-soliton exists only for the defined value 
of a, which can be changed for the fixed absorber properties by 
the variation of A, i. e. by variation of mode cross-section in 
active medium and in absorber or by variation of the absorber 
mirror reflectivity 



Note,that the pulse duration is defined by the formula, which is similar to one 
for Kerr-lens mode locking (see part VII). Let find the pulse duration as the 
function of pump. 

>Energy=2*a0 2 *tp:#energy of quasi-soliton 

fl4 := 
numer(simplify(subs(tp=l/sqrt (gam-alpha), 

subs(Energy=2*a0 2 *tp,fl2)))): 
alpha_sol2 := solve(fl4=0, alpha) :#saturated gain 

fig2 := plot( 
Re(evalf(subs(lambda=0.5,subs( 

{alphamx=0.1,Tr=300,tau=6.25e-4/lambda 2 ,gam=0.01 
},subs(alpha=alpha_sol2[2], 2. 5/sqrt (gam- alpha)))))), 
Pump=0. 0005. .0. 005, axes=boxed,labels=['Pump, a.u.', 'tp,fs'], 
title='Pulse duration versus pump', color=blue): 

display (fig, fig2,view=5.. 20); 
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tp, fs 



Pulse duration versus pump 




0.001 



0.002 „ 0.003 
Pump, a.u. 



0.004 



0.005 



Thus, the Kerr-lensing (lower curve) allows to reduce the pulse duration and to 
generate the sub-10 fs quasi-soliton. 



So, we had demonstrated, that an joint action of the lasing fac- 
tors and coherent absorber objectives to the generation of the 
coherent soliton. The ultrashort pulse in this case has 2 7r area, 
but it is not sec/j-shape pulse (quasi-soliton). The obtained val- 
ues of the pulse duration are placed within interval of 8 - 30 fs. 
The contribution of the self-focusing changes the pulse shape es- 
sentially. In this case, there exists the stable secft-shape quasi- 
soliton with duration, which depends on the absorber mirror 
reflectivity (or ratio of the laser mode cross-section in the active 
medium and absorber). The obtained result is very attractive 
for the elaboration of compact, all-solid-state, "hand-free" fem- 
tosecond lasers 



12 Conclusion 



The powerful computation abilities of Maple 6 allowed to demonstrate the basic 
conceptions of the modern femtosecond technology. The application of these 
conceptions to the Kerr-lens mode locking and mode locking due to coherent 
semiconductor absorber leaded to the new scientific results (see, for example, 
[fl5|), which are very useful for elaboration of the high-stable generators of sub- 
10 fs laser pulses. We can see, that the combination of the symbolical, numerical 
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approaches and programming opens a door for a new opinion on the ultrashort 
pulse generation. This opinion is based on the search of the soliton and quasi- 
soliton states of nonlinear dynamical equations and on the analysis of their 
evolution as the evolution of the breezer's type. 
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